********************************************************************************
**
**  Original code made by Tobias Ferber and Martin Giese (Mandel-92)
**
**  Rewritten and speeded up by Dino Papararo            01-Mar-1997
**
**  FUNCTION
**
**    MandInt -- perform Z = Z^2 + C iteration.
**
**  SYNOPSIS
**
**    WORD MandINT (WORD Iterations,LONG Cre,LONG Cim)
**
**
**  DESCRIPTION
**
**    Work with fixed-point real numbers, x between -8 and 8 (aprox.)
**    being represented by x*2^(28) in a long.
**
**
**    MandInt performs an iteration given by
**
**                             2    2
**    x = r    y = i    x   = x  - y  + r    y   = 2 x  y + i
**     0        0        j+1   j    j         j+1     j  j
**
**
**    2 x y  is calculated via
**
**
**                   2    2    2
**    2 x y = (x + y)  - x  - y
**
**
**    as the squaring of x+y requires only three more multiplicaitons vs.
**    four for a normal multiplication of longs.
**
**    the iterations is performed up to n times, but may stop before that
**    if x^2+y^2 exceeds 2.00
**
**    The value returned is the number of iterations actually done.
**
********************************************************************************

***** d0:Iterations d1:zr d2:zi d3:Quad d4:zr2 d5:zi2 d6:Const d7:Const ********

******************************** a0:Cre a1:Cim *********************************

********************************************************************************


        XDEF  _MandINT

_MandINT:

        movem.l d3-d7,-(sp)     ; save regs


        move.l  d1,a0           ; a0 = Cre
        move.l  d2,a1           ; a1 = Cim
        moveq.l #12,d6          ; d6 = 12
        moveq.l #11,d7          ; d7 = 11


Iter:
        move.l  d1,d3           ; save zr in d3 & calculate zr^2
        bpl.s   Positiv_1
        neg.l   d1

Positiv_1:

        move.w  d1,d4           ; d4   = Rl
        mulu.w  d4,d4           ; d4   = Rl * Rl
        clr.w   d4              ; d4.w = 0
        swap.w  d4              ; d4   = d4 >> 16
        lsr.w   d6,d4           ; d4   = d4 >> 12

        move.w  d1,d5           ; d5   = Rl
        swap.w  d1              ; d1.w = Rh
        mulu.w  d1,d5           ; d5   = Rh * Rl
        lsr.l   d7,d5           ; d5   = d7 >> 11

        add.l   d5,d4           ; d4   = ((Rl * Rl) >> 28) + ((2 * Rl * Rh) >> 12)

        mulu.w  d1,d1           ; d1   = Rh * Rh
        asl.l   #4,d1           ; d1   = d4 << 4
        bvs.s   Exit            ; Check Overflow

        add.l   d1,d4           ; d4 = d4 + ((Rh * Rh) << 4)
                                ; zr * zr


        add.l   d2,d3           ; add zi in d3 & calculate zi^2



        tst.l   d2
        bpl.s   Positiv_2
        neg.l   d2

Positiv_2:

        move.w  d2,d5
        mulu.w  d5,d5
        clr.w   d5
        swap.w  d5
        lsr.w   d6,d5

        move.w  d2,d1
        swap.w  d2
        mulu.w  d2,d1
        lsr.l   d7,d1

        add.l   d1,d5

        mulu.w  d2,d2
        asl.l   #4,d2
        bvs.s   Exit

        add.l   d2,d5           ; d5 = zi * zi



        move.l  d4,d1
        add.l   d5,d1           ; zr2 + zi2
        bvs.s   Exit            ; if overflow exit
        cmpi.l  #$40000000,d1   ; compare (zr2 + zi2) with (4 * 2^28)
        bgt.s   Exit            ; if greater exit



        tst.l   d3              ; calculate (zi + zr)^2
        bpl.s   Positiv_3
        neg.l   d3

Positiv_3:

        move.w  d3,d2
        mulu.w  d2,d2
        clr.w   d2
        swap.w  d2
        lsr.w   d6,d2

        move.w  d3,d1
        swap.w  d3
        mulu.w  d3,d1
        lsr.l   d7,d1

        add.l   d1,d2

        mulu.w  d3,d3
        asl.l   #4,d3
        bvs.s   Exit
        add.l   d3,d2           ; d5 = (zi + zr)^2



        move.l  a0,d1           ; d1  = Cre
        move.l  a1,d3           ; d3  = Cim
        sub.l   d4,d2           ; d2  = (zr + zi)^2 - zr2
        sub.l   d5,d1           ; d1 -= zi2
        add.l   d3,d2           ; d2 += Cim
        add.l   d4,d1           ; d1 += zr2
        sub.l   d5,d2           ; d2 -= zi2

        dbra.w  d0,Iter         ; if Iterations != 0) go to Iter
        moveq.l #0,d0           ; else Iterations = 0

Exit:
        movem.l (sp)+,d3-d7     ; restore regs

        rts                     ; return Iterations

        end                     ; end.