/*
 * MandelVroom 2.0
 *
 * (c) Copyright 1987,1989  Kevin L. Clague, San Jose, CA
 *
 * All rights reserved.
 *
 * Permission is hereby granted to distribute this program's source
 * executable, and documentation for non-comercial purposes, so long as the
 * copyright notices are not removed from the sources, executable or
 * documentation.  This program may not be distributed for a profit without
 * the express written consent of the author Kevin L. Clague.
 *
 * This program is not in the public domain.
 *
 * Fred Fish is expressly granted permission to distribute this program's
 * source and executable as part of the "Fred Fish freely redistributable
 * Amiga software library."
 *
 * Permission is expressly granted for this program and it's source to be
 * distributed as part of the Amicus Amiga software disks, and the
 * First Amiga User Group's Hot Mix disks.
 *
 * contents: this file contains the functions to calculate Mandelbrot and
 * Julia pictures using a special and fast fixed point (scaled ints) format.
 * It has generators for 68000 and 68020 in assembly.
 */

#include "mandp.h"

#include "parms.h"

#define BITS2SHIFT 27

/*
 * 32 bit fixed point generator
 */
MandelbrotInt32( Pict )
  register struct Picture *Pict;
{
  register LONG   i, j, k;
  register SHORT *CountPtr;

  double gap;

  UBYTE ConvFlag;
  LONG  gapx, gapy, startx;

  struct IntPotParms Parms;
  struct RastPort   *Rp;

  if (Pict->Flags & NO_RAM_GENERATE)
    CountPtr = Pict->Counts;
  else
    CountPtr = Pict->Counts + Pict->CurLine*Pict->CountX;

  /* figure out horizontal and verticle distances between points in plane.
   * convert them to fixed point format                                    */

  gapy = (int) (Pict->ImagGap*((double)(1<<BITS2SHIFT)));
  gapx = (int) (Pict->RealGap*((double)(1<<BITS2SHIFT)));

  startx = (int)(Pict->RealLow*((double)(1<<BITS2SHIFT)));

  /*
   * for each point in the image, calculate Mandelbrot
   */
  Parms.C_Imag = (int)(Pict->ImagLow*((double)(1<<BITS2SHIFT)));

  Parms.C_Imag += Pict->CurLine*gapy;

  Parms.ScreenReal = 0;
  Parms.ScreenImag = 0;

  Parms.MaxIteration = Pict->MaxIteration;

  for (i = Pict->CurLine; i < Pict->CountY; i++) {

    Parms.C_Real = startx;

    ConvFlag = 0;

    if ( Pict->Flags & NO_RAM_GENERATE )
      CountPtr = Pict->Counts;

    for (j = 0; j < Pict->CountX; j++) {

      if (*CountPtr == 0) {

        /*
         * if the last pixel was mandelbrot, then use trace table
         */
        if (ConvFlag) {

          k = Int32_Trace_Height( &Parms );
        } else {

          k = Int32_Height( &Parms );
        }

        ConvFlag = k == Pict->MaxIteration;
        *CountPtr = k;
      }
      CountPtr++;

      Parms.C_Real += gapx;

      ChildPause( Pict );
    }
    Parms.C_Imag += gapy;

    CheckEOL( Pict );
  }
} /* Mandelbrot32Int */

#ifdef CHECK_TASK_STACK

int stackerr;
LONG stacklow;
LONG stackhi;
LONG stackovfl;
LONG stackmax;

CheckStack() {
  register struct Task *Task;
  register LONG hi,reg,low;
  register LONG size;

  Task = FindTask(0L);

  hi = Task->tc_SPUpper;
  reg = Task->tc_SPReg;
  low = Task->tc_SPLower;

  size = hi - reg;

  if (reg < low) {

    if (size > stackmax) {
      stackerr  = 1;

      stacklow   = low;
      stackhi    = hi;
      stackovfl  = reg;
      stackmax   = size;
    }
  } else {

    if (stackerr == 0 && size > stackmax) {
      stacklow   = low;
      stackhi    = hi;
      stackovfl  = reg;
      stackmax   = size;
    }
  }
}

PrintStack()
{
  if (stackerr != 0)
    printf("********* STACK ERR ***********\n");

  printf("Low   %08x\nHigh  %08x\nSP     %08x\n",
          stacklow, stackhi, stackovfl );

  printf("MinStack %d\n",stackmax);
}

#endif

/*
 * 32 bit fixed point generator
 */
JuliaInt32( Pict )
  register struct Picture *Pict;
{
  register LONG   i, j, k;
  register SHORT *CountPtr;

  LONG  gapx, gapy, startx;
  UBYTE ConvFlag;

  struct IntPotParms Parms;

  if (Pict->Flags & NO_RAM_GENERATE)
    CountPtr = Pict->Counts;
  else
    CountPtr = Pict->Counts + (Pict->CurLine*Pict->CountX);

  /* figure out horizontal and verticle distances between points in plane.
   * convert them to fixed point format                                    */

  gapx = (int) ( Pict->RealGap * (double) (1 << BITS2SHIFT) );
  gapy = (int) ( Pict->ImagGap * (double) (1 << BITS2SHIFT) );

  /*
   * for each point in the image, calculate Julia
   */
  Parms.ScreenImag  = (int) ( Pict->ImagLow * (double) (1 << BITS2SHIFT) );
  Parms.ScreenImag += Pict->CurLine*gapy;

  startx = (int) ( Pict->RealLow * (double) (1 << BITS2SHIFT) );

  Parms.C_Real = (int) ( Pict->Real * (double) (1 << BITS2SHIFT) );
  Parms.C_Imag = (int) ( Pict->Imag * (double) (1 << BITS2SHIFT) );

  Parms.MaxIteration = Pict->MaxIteration;

  for (i = Pict->CurLine; i < Pict->CountY; i++) {

    Parms.ScreenReal = startx;
    ConvFlag = 0;

    if ( Pict->Flags & NO_RAM_GENERATE )
      CountPtr = Pict->Counts;

    for (j = 0; j < Pict->CountX; j++) {

      if (*CountPtr == 0) {

        /*
         * if the last pixel was mandelbrot, then use trace table
         */

        if ( ConvFlag ) {
          k = Int32_Trace_Height( &Parms );
        } else {
          k = Int32_Height( &Parms );
        }

        ConvFlag = k == Pict->MaxIteration;
        *CountPtr = k;
      }
      CountPtr++;

      Parms.ScreenReal += gapx;

      ChildPause( Pict );
    }
    Parms.ScreenImag += gapy;

    CheckEOL( Pict );
  }
} /* Julia32Int */

/*
 * 32 bit fixed point generator
 */
MandelbrotInt32II( Pict )
  register struct Picture *Pict;
{
  register LONG   i, j, k;
  register struct RastPort *Rp;

  double gap;

  SHORT *CountPtr;
  LONG  gapx, gapy, startx;

  struct IntPotParms Parms;

  if (Pict->Flags & NO_RAM_GENERATE)
    CountPtr = Pict->Counts;
  else
    CountPtr = Pict->Counts + Pict->CurLine*Pict->CountX;

  /* figure out horizontal and verticle distances between points in plane.
   * convert them to fixed point format                                    */

  gapy = (int) (Pict->ImagGap*((double)(1<<BITS2SHIFT)));
  gapx = (int) (Pict->RealGap*((double)(1<<BITS2SHIFT)));

  startx = (int)(Pict->RealLow*((double)(1<<BITS2SHIFT)));

  /*
   * for each point in the image, calculate Mandelbrot
   */
  Parms.C_Imag = (int)(Pict->ImagLow*((double)(1<<BITS2SHIFT)));
  Parms.C_Imag += Pict->CurLine*gapy;

  Parms.ScreenReal = 0;
  Parms.ScreenImag = 0;

  Parms.MaxIteration = Pict->MaxIteration;

  for (i = Pict->CurLine; i < Pict->CountY; i++) {

    Parms.C_Real = startx;

    if ( Pict->Flags & NO_RAM_GENERATE )
      CountPtr = Pict->Counts;

    for (j = 0; j < Pict->CountX; j++) {

      /*
       * if the last pixel was mandelbrot, then use trace table
       */

      k = Int32_Height_Fast( &Parms );

      *CountPtr++ = k;

      Parms.C_Real += gapx;

      ChildPause( Pict );
    }
    Parms.C_Imag += gapy;

    CheckEOL( Pict );
  }
} /* MandelbrotInt32II */

/*
 *  This code does mandelbrot without any trace lookup.
 *  It is the fastest 68000 generator in the house.
 */
Int32_Height( Parms )

  struct IntPotParms *Parms;
{
  LONG Height;

  register struct IntPotParms *P = Parms;

  register LONG cura, curb, cura2, curb2;

#asm
height  equ -4
Bits2Shift equ  5
;
;
;  d1 - BITS2SHIFT
;  d2 - k
;  d4 - a
;  d5 - b
;  d6 - a2
;  d7 - b2
;
screenr  equ 0
screeni  equ 4
curx     equ 8
cury     equ 12
maxi     equ 16
;
   move.l   #Bits2Shift,d1
   move.l   screenr(a2),d4
   move.l   screeni(a2),d5
   move.w   maxi(a2),d2
   bra      Fposa2
;
FKLoop
;
;  cura = cura2 - curb2 + curx;
   exg      d6,d4       ; exchange cura and cura2
   sub.l    d7,d4       ; subtract curb
   add.l    curx(a2),d4 ; add curx
;
;  curb = cura * curb >> 12;
   move.l   d6,d7      ; get copy of op1 sign bit
   bpl      Fpos1       ; get absolute value of op1
   neg.l    d6
Fpos1
   eor.l    d5,d7      ; calculate result sign
   tst.l    d5         ; get absolute value of op2
   bpl      Fpos2
   neg.l    d5
Fpos2
   move.l   d6,d0      ; get a copy of op1
   swap     d0         ; get high half of op1
   move.w   d0,d7      ; save a copy of high half
   mulu     d5,d0      ; multiply op2 low by op1 high
   clr.w    d0         ;  clear least significant part
   swap     d0         ;  put it in it's place
   swap     d5         ; get high half of op2
   mulu     d5,d6      ; multiply op2 high with op1 low
   clr.w    d6         ;  clear least significant part
   swap     d6         ;  put it in its place
   mulu     d7,d5      ; multiply op2 high by op1 high
   add.l    d0,d5      ; add partial results
   add.l    d6,d5      ; add partial results
   tst.l    d7         ; is the result negative?
   bpl      Fpos3
   neg.l    d5         ; yes, better negate it.
Fpos3
   asl.l    d1,d5      ; now, rescale it.
;
;  curb += curb + cury;
   add.l    d5,d5      ; double it and add cury
   add.l    cury(a2),d5
Fposa2
;
;  cura2 = cura * cura;
   move.l   d4,d0      ; get absolute value of a in d0
   bpl      Fposa
   neg.l    d0
Fposa
   move.l   d0,d6      ; copy absolute value into d6
   swap     d6         ; get high part in d6
   mulu     d6,d0      ; multiply high and low destroying low
   clr.w    d0         ; clear the least significant part
   swap     d0         ; put most sig. part in low half
   mulu     d6,d6      ; multiply high and high destroing high
   add.l    d0,d6      ; add in lower half twice
   add.l    d0,d6
   asl.l    d1,d6      ; get radix point back in correct place
   bvs      Fbailout
;
;  curb2 = curb * curb;
   move.l   d5,d0      ; get absolute value of a in d0
   bpl      Fposb
   neg.l    d0
Fposb
   move.l   d0,d7      ; copy absolute value into d7
   swap     d7         ; get high part in d7
   mulu     d7,d0      ; multiply high and low destroying low
   clr.w    d0         ; clear the least significant part
   swap     d0         ; put most sig. part in low half
   mulu     d7,d7      ; multiply high and high destroing high
   add.l    d0,d7      ; add in lower half twice
   add.l    d0,d7
   asl.l    d1,d7      ; get radix point back in correct place
   bvs      Fbailout
;
   move.l   d6,d0      ; if (cura2 + curb2 >= 4) goto bailout;
   add.l    d7,d0
   bvs      Fbailout
;
   dbra     d2,FKLoop
;  addq     #1,d2
Fbailout
   move.w   maxi(a2),d0
   sub.w    d2,d0
   sub.w    #1,d0
   ext.l    d0
#endasm
   ;;
}

Int32_Trace_Height( Parms )

  register struct IntPotParms *Parms;
{
  SHORT  k,TLen;
  LONG  OldA[17];

  register LONG cura, curb, cura2, curb2;

#asm
;
;
;  d1 - BITS2SHIFT
;  d2 - k
;  d4 - a
;  d5 - b
;  d6 - a2
;  d7 - b2
;  a0 - Trace Table pointer
;
Trace equ 16
OldA  equ -132
k     equ -2
TLen  equ -4

   move.l   #Bits2Shift,d1
   move.l   screenr(a2),d4
   move.l   screeni(a2),d5
;
   lea      OldA(a5),a0          ; Set up Trace table pointer
   move.w   #Trace,d2            ; Set up Trace table len
;
   move.w   #-1,k(a5)
   bra      TPosa2               ; branch in to middle to get a2 = a * a
;
TLoop
;
;  cura = cura2 - curb2 + curx;
   exg      d6,d4       ; exchange cura and cura2
   sub.l    d7,d4       ; subtract curb
   add.l    curx(a2),d4 ; add curx
;
;  curb = cura * curb >> 12;
   move.l   d6,d7      ; get copy of op1 sign bit
   bpl      TPos1       ; get absolute value of op1
   neg.l    d6
TPos1
   eor.l    d5,d7      ; calculate result sign
   tst.l    d5         ; get absolute value of op2
   bpl      TPos2
   neg.l    d5
TPos2
   move.l   d6,d0      ; get a copy of op1
   swap     d0         ; get high half of op1
   move.w   d0,d7      ; save a copy of high half
   mulu     d5,d0      ; multiply op2 low by op1 high
   clr.w    d0         ;  clear least significant part
   swap     d0         ;  put it in it's place
   swap     d5         ; get high half of op2
   mulu     d5,d6      ; multiply op2 high with op1 low
   clr.w    d6         ;  clear least significant part
   swap     d6         ;  put it in its place
   mulu     d7,d5      ; multiply op2 high by op1 high
   add.l    d0,d5      ; add partial results
   add.l    d6,d5      ; add partial results
   tst.l    d7         ; is the result negative?
   bpl      TPos3
   neg.l    d5         ; yes, better negate it.
TPos3
   asl.l    d1,d5      ; now, rescale it.
;
;  curb += curb + cury;
   add.l    d5,d5      ; double it and add cury
   add.l    cury(a2),d5
TPosa2
;
;  cura2 = cura * cura;
   move.l   d4,d0      ; get absolute value of a in d0
   bpl      TPosa
   neg.l    d0
TPosa
   move.l   d0,d6      ; copy absolute value into d6
   swap     d6         ; get high part in d6
   mulu     d6,d0      ; multiply high and low destroying low
   clr.w    d0         ; clear the least significant part
   swap     d0         ; put most sig. part in low half
   mulu     d6,d6      ; multiply high and high destroing high
   add.l    d0,d6      ; add in lower half twice
   add.l    d0,d6
   asl.l    d1,d6      ; get radix point back in correct place
   bvs      bailout
;
;  curb2 = curb * curb;
   move.l   d5,d0      ; get absolute value of a in d0
   bpl      TPosb
   neg.l    d0
TPosb
   move.l   d0,d7      ; copy absolute value into d7
   swap     d7         ; get high part in d7
   mulu     d7,d0      ; multiply high and low destroying low
   clr.w    d0         ; clear the least significant part
   swap     d0         ; put most sig. part in low half
   mulu     d7,d7      ; multiply high and high destroing high
   add.l    d0,d7      ; add in lower half twice
   add.l    d0,d7
   asl.l    d1,d7      ; get radix point back in correct place
   bvs      bailout

   move.l   d6,d0      ; if (cura2 + curb2 >= 4) goto bailout;
   add.l    d7,d0
   bvs      bailout

   move     d0,(a0)+   ; save magnitude in trace table

   addq.w   #1,k(a5)   ; are we out of iterations?
   move.w   maxi(a2),d0
   cmp.w    k(a5),d0
   ble      GotIt
   dbra     d2,TLoop   ; nope, so try again till trace table full

   move.w   -(a0),d0   ; get last entry in the trace

   move.w   #Trace-1,d2  ; set up length and address for comparison loop
   lea      OldA(a5),a0

TCLoop
   cmp      (a0)+,d0
   beq      GotIt       ; did we find it?
   dbra     d2,TCLoop

   lea      OldA(a5),a0 ; no match, so prepare for the next round
   move     d0,(a0)+    ; move last entry of the table

   move.w   #Trace,d2
   bra      TLoop

GotIt
   move.w   maxi(a2),k(a5)
bailout
   move.w   k(a5),d0
   ext.l    d0
#endasm
   ;;
}

/*
 * 32 bit fixed point generator
 */
Int32_Height_Fast( Parms )
  struct IntPotParms *Parms;
{
  register struct IntPotParms *P = Parms;
  register LONG cura,curb,cura2,curb2;
  /*
   *  This is the fastest generator in the house.
   */
#asm
   machine mc68020
fscreenr  equ 0
fscreeni  equ 4
fcurx     equ 8
fcury     equ 12
fmaxi     equ 16

Zcurx      equ   8
Zcury      equ  12
;
;  d1 - BITS2SHIFT
;  d2 - k
;  d4 - a
;  d5 - b
;  d6 - a2
;  d7 - b2
;
;  cura = curb = curc = curd = 0;
   move.w   maxi(a2),d1
   move.l   #27,d2
   move.l   fscreenr(a2),d4
   move.l   fscreeni(a2),d5
   move.l   fcurx(a2),a0
   move.l   fcury(a2),a1
   bra      Zpos2

ZLoop
;
;  cura = cura2 - curb2 + curx;
;
   exg      d6,d4
   sub.l    d7,d4
;   move.l   a0,d0
   add.l    a0,d4
;
;  curb = cura * curb
;
   muls.l   d6,d0:d5
   asl.l    #5,d0
   lsr.l    d2,d5
   or.l     d0,d5
;
;  curb += curb + cury;
;
   add.l    d5,d5
;   move.l   a1,d0
   add.l    a1,d5
Zpos2
;
;  cura2 = cura * cura;
;
   move.l   d4,d6
   muls.l   d6,d0:d6
   asl.l    #5,d0
   bvs      Zbailout
   lsr.l    d2,d6
   or.l     d0,d6
;
;  curb2 = curb * curb;
;
   move.l   d5,d7
   muls.l   d7,d0:d7
   asl.l    #5,d0
   bvs      Zbailout
   lsr.l    d2,d7
   or.l     d0,d7
;
   move.l   d6,d0
   add.l    d7,d0
   bvs      Zbailout
;
;  cmp.l    #536870912,d0
;   bge      Zbailout
;
   dbra     d1,ZLoop
;
;  addq     #1,d1
Zbailout
   move.w   fmaxi(a2),d0
   sub.w    d1,d0
   sub.w    #1,d0
   ext.l    d0
#endasm
   ;
}