********************************************************************
*              Routine Mandelbrot/Julia pour HiSoft-Basic
* Emploi:
*       CALL LOC Routine&,VARPTR(Args%(0)),VARPTR(ColorTable%(0))
*       avec Arg%()        = tableau des arguments d'entrée
*            ColorTable%() = table des couleurs
*            (voir à la fin du programme)
* Si le calcul est interrompu, la routine place un 'X' en tête de Args%(0)
*********************************************************************

* Modifications sur la phase 2.0 :
*    l'argument d'entrée RASTPORT est remplacé par WINDOW 
*    l'arrêt n'est déclenché que si l'écran est en "FirstScreen"
*    débogage du remplissage de Tradu (avr 91)
*    ajout d'un pointeur spécial

*   OPT O+

ExecBase=4
OpenLibrary  =-408
CloseLibrary =-414
AllocMem     =-198
FreeMem      =-210
AvailMem     =-216

ReadPixel    =-318
WritePixel   =-324
SetAPen      =-342
SetDrMd      =-354

Open         =-30
Close        =-36
Read         =-42
Write        =-48
Seek         =-66
IoErr        =-132

PUSH:   macro
       lea \3(pc),a0
       move.\1 \2,(a0)
       ENDM

***************** en tête pour mise au point ********************** 

* (Cet en-tête devra être supprimé avant la compilation pour le Basic)

debut:
       lea    dosname(pc),a1
       move.l 4,a6
       jsr    -408(a6)
       PUSH   l,d0,dosbase
       lea    gfxname(pc),a1
       move.l 4,a6
       jsr    -408(a6)
       PUSH   l,d0,gfxbase
       lea    intuiname(pc),a1
       jsr    -408(a6)
       PUSH   l,d0,intuibase
       move.l d0,a6
       lea    NScreen(pc),a0
       jsr    -198(a6)
       PUSH   l,d0,myscreen1
       add.l  #$54,d0
       PUSH   l,d0,myrp
       lea    Nwindow(pc),a0
       jsr    -204(a6)             OpenWindow(a6)
       PUSH   l,d0,mywindow1
       PUSH   l,d0,mywindow2

       move.l myscreen1(pc),d0
       add.l  #44,d0
       move.l d0,a0
       lea    mycolors(pc),a1
       move.l #32,d0
       move.l gfxbase(pc),a6
       jsr    -192(a6)            LoadRGB4

       lea    ColorTable(pc),a0
       move.l a0,-(sp)
       lea    Arguments(pc),a0
       move.l a0,-(sp)
       bsr    start
       movem.l (sp)+,d0/d1

wait:  btst #6,$BFE001
       bne.s  wait

       move.l intuibase(pc),a6
       move.l mywindow1(pc),a0
       jsr    -72(a6)               CloseWindow(a6)

       move.l myscreen1(pc),a0
       jsr    -66(a6)
       move.l 4,a6
       move.l intuibase(pc),a1
       jsr    -414(a6)
       move.l gfxbase(pc),a1
       jsr    -414(a6)
       rts

Large=320
Haut=200
ViewMode=$0000

NScreen:    dc.w 0,0,Large,Haut,5,$0102,ViewMode,15,0,0,0,0,0,0
Nwindow:    dc.w 0,0,Large,Haut,$102
            dc.l 0,1,0,0,0  
myscreen1:  dc.l 0,0,0,0
            dc.w 15 

myrp:       dc.l 0
mywindow1:  dc.l 0
mycolors:   dc.w $0000,$0999,$0444,$0555,$0666,$0777,$0888,$0999
            dc.w $0AAA,$0BBB,$0CCC,$0DDD,$0EEE,$0FFF,$0FEE,$0FDD
            dc.w $0FCC,$0FBB,$0FAA,$0F99,$0F88,$0F77,$0F66,$0F55
            dc.w $0E64,$0D73,$0C82,$0B91,$0AA0,$0BB0,$0CC0,$0DD0

* Arguments d'entrée
Max=60


Arguments:  dc.b 'P',32,1,1        Type,NbBits, ModeCalcul,Iter
*       dc.l -2147483648,2147483333          X1
*       dc.l  254898349,254898666           Y1
*       dc.l -1451732491,0          X2
*       dc.l -266915032,0           Y2
       dc.l $80000000,0,$4ccc0000,0       X1,Y1  pour l'ensemble général
       dc.l $4ccc0000,0,$b3330000,0       X2,Y2
*    dc.l $a51eb85b,0,$0866666c,0,$b1eb852b,0,$fecccccc,0
*       dc.w  160,128                NX,NY
       dc.w  Large,Haut            NX,NY
       dc.w  185,70,258            XR1,YR1,XR2
       dc.w  Large,Haut            NX0,NY0
       dc.l  0,0,0                 DX(2),DY
*       dc.l  -205000000,0                  XS
       dc.l  0,0        -405000000,0                  XS
       dc.l  0,0         670000000,0                  YS
mywindow2       dc.l 0                    Window
       dc.w Max,30               NMax,NbCouleurs
       dc.l 0                    => table des couleurs (NbCouleurs entrées)
       dc.l  IterFile 
       dcb.w 3,0                espace de travail (contient aussi les 2 mots 
       dcb.w 9,0                                   de IterFIle)
       dc.l $80000000,0            X1J
       dc.l $60000000,0            Y1J 
       dc.l $7fffffff,0            X2J
       dc.l $a0000000,0            Y2J


ColorTable: dc.w  0      0/1 pour couleurs normales/ditherées
            dc.w  2,0, 3,1, 4,2, 5,3, 6,4, 7,5, 8,6, 9,7
*            dc.w  28,0, 27,1, 26,2, 25,3, 24,4, 23,5, 22,6, 21,7
            dc.w  10,8, 11,9, 12,10, 13,11, 14,12, 15,13, 16,14, 17,15
            dc.w  18,16,19,17,20,18, 21,19, 22,20, 23,21, 24,22, 25,23
            dc.w  26,24,27,25,28,26, 29,27, 30,28, 31,Max-1 

IterFile:  dc.b 'df1:Iter',0
      even

************************* fin en-tête *****************************



************************  ARGUMENTS  *******************************

*    Type="M" (pour Mandelbrot), "J" (Julia), "1" à "5" (Mandelbrot incomplets)
*         "r" (simple recoloration)
*         "R" (recoloration avec redistribution des couleurs)    (1octet)
*
*    NbBits         : 32,48 ou 64 bits pour le calcul            (1octet)
*
*    ModeCalcul = 0 : agrandissement normal du rectangle défini par
*                     XR1,YR1,XR2 dans l'image de base X1,Y1,X2,Y2,NX,NY
*                     quel que soit le type
*                     (X1,Y1,X2,Y2 sont redéfinis)
*               = 1 : calcul de l'image défini par X1,Y1,X2,Y2,NX,NY 
*                     en Mandelbrot pur, ou bien, en Julia ou Mandelbrot
*                     avec source, calcul de la source XS,YS à partir   
*                     de X1,X2,Y1,Y2 et XR1,YR1, et redéfinition de X1,..
*                     pour une image plein format.               (1octet)
*
*    Iter = 0 :  pas de fichier "iter" et pas de recoloration
*         = 1 :  fichier "iter" en ram:
*         = 2 :  fichier "iter" transmis à l'appel à l'offset Work
*                    dans  le tableau Args()                      (1octet)                  
*
* X1,Y1,X2,Y2 : coordonnées sur 64 bits du cadre de l'image de référence
*               (image à agrandir, ou bien image Mandelbrot de référence 
*                dans le cas d'une nouvelle image de type Julia)  (8 longs)
*               
* NX, NY      : nombre de pixels horizontaux et verticaux         (2 mots)
*               de l'image à construire
*
* XR1,YR1,XR2 : en agrandissement, coordonnées (pixels) du cadre à
*               agrandir. La 4ème coordonnée est fixée avec l'hypothèse
*               que l'image de départ a un rapport de dimension NX0/NY0
*               et que l'on conserve ce rapport dans l'agrandissement.
*               En mode Julia nouveau, XR1 et YR1 sont les coordonnées
*               du point source, à prendre dans le cadre (X1,Y1)-(X2,Y2)
*               de référence.                                      (3 mots)
*
* NX0,NY0     : dimensions (pixels) de l'image de référence. EFFACES EN
*               SORTIE                                             (2 mots)
* 
* DX2,DY      : zone de calcul interne                             (3 longs) 
* 
* XS,YS       : coordonnées (64 bits) du point source              (4 longs)
*
* Window      : adresse de la structure Window à remplir (donnée par 
*               la fonction WINDOW(7) en Basic)                    (1 long) 
*
* NMax        : nombre maximum d'itérations                        (1 mot)
*
* NbCol       : nombre de couleurs utilisées dans la palette (numéros 1
*               à NbCol)                                           (1 mot) 
*
* ColTab      : champ interne (reprise de l'adresse de la table de translation)
*                                                                  (1 long)
* Work et Work2 : 28 octets utilisés comme zone de travail interne.
*               En entrée, lorsque 'Iter'=2, on transmet l'adresse du nom du
*               fichier Iter à l'adresse Work.
*               En sortie, on trouve NMin à cette adresse
*
* X1J,Y1J,X2J,Y2J : cadre (en 64 bits) des images de type Julia, soit en
*               entrée (nouvelle source) soit en sortie (agrandissement).
*               Dans ce dernier cas, copie des X1...Y2             (8 longs)

 
* La table des couleurs se compose d'un premier mot indiquant si les 
* couleurs sont normales ou dithérées. Ensuite :

*   Si les couleurs sont normales, il y a ensuite NbCouleurs couples de
*   mots   [registre couleur, max d'itérations pour cette couleur] ; il
*   peut y avoir des répétitions de couleur.

*   Si les couleurs sont dithérées (premier mot non nul), philosophie
*   non encore fixée

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





*
* a5 va pointer sur la liste des arguments d'entrée
* Vérifier que le pointeur de fenêtre existe (sinon, sortir)

start:
       movem.l d0-d7/a0-a6,-(sp)
       move.l 64(sp),a5          
       move.l 68(sp),ColTab(a5)
       move.l Window(a5),d0
       beq.s  ouf

* On va en déduire les adresses d'écran et de rastport
       move.l d0,a1
       lea    myscreen(pc),a0
       move.l $2E(a1),(a0)+        screen
       move.l $32(a1),(a0)         rastport

       lea    dosname(pc),a1
       move.l 4,a6
       jsr    OpenLibrary(a6)
       PUSH   l,d0,dosbase
       beq.s  ouf
       lea    gfxname(pc),a1
       jsr    OpenLibrary(a6)
       PUSH   l,d0,gfxbase
       beq.s  exdos
       lea    intuiname(pc),a1
       jsr    OpenLibrary(a6)
       PUSH   l,d0,intuibase
       beq.s  exgfx
aaa       move.l d0,a6

       bsr    trace

bbb
* On ferme      

       move.l 4,a6
       move.l intuibase(pc),a1
       jsr    CloseLibrary(a6)
exgfx: move.l gfxbase(pc),a1
       jsr    CloseLibrary(a6)
exdos: move.l dosbase(pc),a1
       jsr    CloseLibrary(a6)
ouf:   
       movem.l (sp)+,d0-d7/a0-a6
       rts

intuiname:  dc.b 'intuition.library',0
gfxname:    dc.b 'graphics.library',0
dosname:    dc.b 'dos.library',0
  even
intuibase:  dc.l 0
gfxbase:    dc.l 0
dosbase:    dc.l 0
myscreen:   dc.l 0
RastPort:   dc.l 0


* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


* Calcul d'une image de Mandelbrot/Julia
* a5 pointe sur la liste des arguments d'entrée
* a4 va pointer sur les arguments internes de la routine de calcul
* a3 va pointer sur un tampon "Ligne"
* Coordonnées d'écran = d6,d7

trace: 

* ouverture des buffers Histo, Tradu, Ligne et du fichier 'iter'
       bsr    OpenAll
       beq    fin
* initialisation de PixTable
       bsr    InitPix
* si Args(0) est 'r' ou 'R', on ne fait qu'une recoloration
       move.b (a5),d0
       and.b  #$DF,d0
       cmp.b  #'R',d0
       beq    recolore      



* Calcul de X1,X2,...DX,DY selon le mode de calcul
* Si ModeCalcul=1 et mode<>"M", redéfinition du point source 
* Si ModeCalcul=0, on calcule un agrandissement
* On utilise un espace de travail de 14 mots (en a5+Work)

       lea    Xcur0(pc),a4

       tst.b  ModeCalcul(a5)
       beq.s  NewXY
       cmp.b  #'M',(a5)
       beq    CalculDXDY

* X2 = X1 + (X2-X1)*XR2/NX0 et X1 = X1 + (X2-X1)*XR1/NX0 
* On met (X2-X1) à l'adresse Work+0, sur 64 bits 

NewXY:
       move.l a5,a0
       move.l a5,a1
       adda.l #Work+8,a1 
       adda.l #X1+8,a0                a0 juste aprés X1
       move.l X2+4(a5),Work+4(a5)
       move.l X2(a5),Work(a5)
       move   #0,ccr
       subx.l -(a0),-(a1) 
       subx.l -(a0),-(a1)             a1 maintenant sur Work+0
       move.l a5,a2
       move.l a5,a3
       adda.l #X1,a2
       adda.l #X2,a3
       move.w XR2(a5),d1
       move.w NX0(a5),d0
       bsr    MulDivAdd               recalcule X2               

       move.l a5,a1
       move.l a5,a2
       adda.l #X1,a2
       move.l a2,a3
       adda.l #Work,a1
       move.w XR1(a5),d1
       move.w NX0(a5),d0
       bsr    MulDivAdd               recalcule X1

* Y1 = Y2 + (Y1-Y2)*(NY0-YR1)/NY0 et Y2 = Y1 + (Y1-Y2)*(XR2-XR1)/NX0 
* On met (Y2-Y1) à l'adresse Work+0, sur 64 bits 
       move.l a5,a0
       move.l a5,a1
       adda.l #Work+8,a1 
       adda.l #Y2+8,a0                a0 juste aprés Y1
       move.l Y1+4(a5),Work+4(a5)
       move.l Y1(a5),Work(a5)
       move   #0,ccr
       subx.l -(a0),-(a1) 
       subx.l -(a0),-(a1)             a1 maintenant sur Work+0
       move.l a5,a2
       move.l a5,a3
       adda.l #Y2,a2
       adda.l #Y1,a3
       move.w NY0(a5),d1
       sub.w  YR1(a5),d1
       move.w NY0(a5),d0
       bsr    MulDivAdd               recalcule Y2 

       move.l a5,a1
       move.l a5,a2
       move.l a5,a3
       adda.l #Y1,a2
       adda.l #Y2,a3
       adda.l #Work,a1
       move.w XR2(a5),d1
       sub.w  XR1(a5),d1
       move.w NX0(a5),d0
       bsr    MulDivSub               recalcule Y1

*  En Julia, si ModeCalcul=1, alors  XS=X1 et YS=Y1, et redéfinition 
*     du cadre de l'image d'après X1J...Y2J
*  En Julia, si ModeCalcul=0 (agrandissement), recopie du cadre X1...
*     en X1J...
*  

       move.l a5,a0
       move.l a5,a1
       adda.l #X1,a0
       cmp.b  #'M',(a5)
       beq.s  CalculDXDY
       tst.b  ModeCalcul(a5)
       beq.s  copieXY
       adda.l #XS,a1
       moveq  #3,d0
nxy1   move.l (a0)+,(a1)+            transfert X1,Y1 -> XS,YS
       dbf    d0,nxy1
       adda.l #Y2J-YS,a1             a1->juste après Y2J
       adda.l #16,a0                 a0->juste après Y2 
       moveq  #7,d0
nxy2   move.l -(a1),-(a0)            transfert X1J...Y2J -> X1...Y2
       dbf    d0,nxy2
       bra.s  CalculDXDY

copieXY:
       adda.l #X1J,a1
       moveq  #7,d0
nxy3   move.l (a0)+,(a1)+            transfert X1... -> XIJ...
       dbf    d0,nxy3 


* On va agrandir le rectangle (X1,Y1)-(X2,Y2). On calcule DX et DY (64 bits)
* en divisant par NX et NY (image cible)    

CalculDXDY:
       move.l a5,a0
       move.l a5,a1
       adda.l #Work+6,a0
       adda.l #X1+4,a1
       move.w #0,Work(a5)
       move.l X2+4(a5),Work+6(a5)
       move.l X2(a5),Work+2(a5)
       move.l (a1),d0
       sub.l  d0,(a0)
       subx.l -(a1),-(a0)
       move.w #0,-(a0)
       move.l a5,a1
       adda.l #DX,a1
       move.w NX(a5),d0
       bsr    Div65

       move.l a5,a0
       move.l a5,a1
       adda.l #Work+6,a0
       adda.l #Y2+4,a1
       move.w #0,Work(a5)
       move.l Y1+4(a5),Work+6(a5)
       move.l Y1(a5),Work+2(a5)
       move.l (a1),d0
       sub.l  d0,(a0)
       subx.l -(a1),-(a0)
       move.w #0,-(a0)
       move.l a5,a1
       adda.l #DY,a1
       move.w NY(a5),d0
       bsr    Div65




* Ecriture de l'en-tête du fichier "Iter" : 'ITER',NX,NY. On se sert du
*  buffer 'Ligne'  comme tampon

       move.l IterHandle(pc),d1
       beq.s  debutImage
       move.l Ligne(pc),a0
       move.l a0,d2
       move.l #'ITER',(a0)+
       move.w NX(a5),(a0)+
       move.w NY(a5),(a0)
       move.l #8,d3
       jsr    Write(a6)
       tst.l  d0
       bmi    SortieSec             Write error!

* Calcul de l'image 
* Les résultats de calcul sont stockés dans le buffer 'Ligne' jusqu'à
*  ce qu'on y ait mis 1 ou 32 lignes (selon que Iter<2 ou non).
* On se sert de 'CompteLignes(a5)' comme compteur de lignes
 
       move.l #0,SetFlag(a5)        Ce flag sera mis à 1 dans le Set principal

debutImage
       move.l Y1(a5),8(a4)          chargement Y1
       move.l Y1+4(a5),12(a4)
       moveq  #0,d7                 ligne écran
       move.w NbLinBuffer(pc),CompteLignes(a5)
       move.l Ligne(pc),a3          initialisation du tampon de ligne

NewLine:
       move.l X1(a5),(a4)           chargement X1
       move.l X1+4(a5),4(a4)
       moveq  #0,d6                 numéro pixel horizondal

NewPixel:
       movem.l d6/d7/a3,-(sp)
       bsr     Calcul           le nombre de tours revient en d0
       movem.l (sp)+,d6/d7/a3
       btst    #6,$BFE001       stopper si on clique
       bne.s   insc
       move.l  intuibase(pc),a6
       move.l  $3C(a6),d1         le FirstScreen du moment
       cmp.l   myscreen(pc),d1    On ne s'arrête que si ce FirstScreen
       bne.s   insc               est bien notre écran

* On demande confirmation
       movem.l d0-d7/a0-a6,-(sp)
       lea     intuitxt0(pc),a1
       move.l  a1,a2
       adda.l  #20,a2              a2 pointe sur intuitxt1
       move.l  a2,a3
       adda.l  #20,a3              et a3 sur intuitxt2
       lea     text0(pc),a0
       move.l  a0,12(a1)           linkages des textes
       lea     text1(pc),a0
       move.l  a0,12(a2)
       lea     text2(pc),a0
       move.l  a0,12(a3)
       move.l  Window(a5),a0
       move.l  intuibase(pc),a6
       moveq   #$20,d0                 GADGETDOWN
       move.l  d0,d1                 GADGETDOWN
       move.l  #170,d2               largeur
       moveq   #60,d3                hauteur
       jsr     -348(a6)              AutoRequest
       tst.l   d0
       movem.l (sp)+,d0-d7/a0-a6
       beq.s   insc
SortieSec:
       move.b  #'X',(a5)     
       bra     fin


* Inscription
insc:  move.l Histo(pc),a0
       move.w d0,d1
       add.w  d1,d1                  
       add.w  #1,0(a0,d1.w)
       move.w d0,(a3)+          inscription (mot) dans "Ligne"

       move.l Tradu(pc),a0
       move.w 0(a0,d1.w),d0
       bsr    WritePix

       addq.w  #1,d6
       cmp.w   NX(a5),d6
       beq.s   FinLigne
       move.l  (a4),d0
       move.l  4(a4),d1
       move.l  DX(a5),d2
       add.l   DX+4(a5),d1
       addx.l  d2,d0
       move.l  d0,(a4)
       move.l  d1,4(a4)  
       bra     NewPixel
FinLigne:
       move.l  8(a4),d0  
       move.l  12(a4),d1
       move.l  DY(a5),d2
       sub.l   DY+4(a5),d1
       subx.l  d2,d0
       move.l  d0,8(a4)
       move.l  d1,12(a4)  
       addq.w  #1,d7
       cmp.w   NY(a5),d7            fin de l'image?
       beq.s   transfert            si oui, dernier transfert
       sub.w   #1,CompteLignes(a5)
       bne.s   FinLigne1

* recopie du tampon Ligne dans Iter
transfert:
       move.l  IterHandle(pc),d1
       beq.s   transf1              sauf s'il n'y a pas de Iter...
       move.l  Ligne(pc),d2
       move.l  a3,d3
       sub.l   Ligne(pc),d3
       move.l  dosbase(pc),a6
       jsr     Write(a6)
       tst.l   d0                   
       bmi     SortieSec            write error! 
transf1:
       move.l  Ligne(pc),a3         réinitialisation du tampon 'Ligne'
       move.w  NbLinBuffer(pc),CompteLignes(a5)

FinLigne1:
       cmp.w   NY(a5),d7            fin de l'image?
       bne     NewLine

recolore:
* sauf s'il n'y a pas de ram:iter
       move.l  IterHandle(pc),d1
       beq     fin  

* La recoloration dépend de (a5) (le 1er octet du tableau Args) 
* Si (a5)='r' on suit la table Tradu initiale. On saute à 'Colorie'
* avec d5<>0 
* Si (a5)='M','J' ou 'R', on reforme la table Tradu à partir de
* l'histogramme. Celui-ci doit être recréé si (a5)='R'. Cela s'obtient
* en sautant à 'Colorie' avec d5=0. 

       moveq  #0,d5
       cmp.b  #'R',(a5)
       beq    Colorie
       moveq  #1,d5
       cmp.b  #'r',(a5)
       beq    Colorie 


* Remplissage de la table de conversion et de la table des couleurs
* a0 -> Histo     a6-> Hist(Nmax)  (a0 ne devra pas atteindre a6)
* a1 -> Table des couleurs. a1 pointera sur les nombres ; les couleurs
*       correspondantes seront -2(a1)
* a2 -> Tradu
* d5 : nombre de couleurs restant à répartir
* d7 : nombre de points restant à colorier


rc0:   move.l  Histo(pc),a0
       move.l  ColTab(a5),a1
       adda.l  #4,a1
       move.l  Tradu(pc),a2
       move.w  NbCol(a5),d5     d5=Nombre de couleurs restant à répartir

* Calcul du nombre initial de points hors NMax, dans d7
       moveq   #0,d7
       moveq   #0,d2
       move.w  NMax(a5),d0
       subq.w  #1,d0
       move.l  a0,a6
rc1:   move.w  (a6)+,d2
       add.l   d2,d7
       dbf     d0,rc1           a6 pointe maintenant vers Histo(NMax)
       moveq   #-1,d6           compteur de tours d'iterations   

rc2:   move.l  d7,d0
       divu    d5,d0            piège d'overflow???
       and.l   #$FFFF,d0
       moveq   #-1,d3
* d0 contient le nombre moyen de points par couleur restante; d3 va servir
* à détecter la 1ère itération  

rc3    tst.l   d7
       beq.s   rcEX
       cmp.l   a0,a6
       beq.s   rcEX
       move.w  d0,d1           mémorisation
       move.w  (a0)+,d2
       addq.w  #1,d3
       addq.w  #1,d6
       sub.l   d2,d7           remise à jour de NTOTAL
       move.w  -2(a1),(a2)+    remplissage TRADU
       sub.l   d2,d0   
       bpl.s   rc3
* Si d3=0 (1ère itération) poursuivre
       tst.w   d3
       beq.s   rc4
* Sinon, comparer d1 et abs(d0)
       neg.l   d0
       cmp.w   d0,d1
       bcc.s   rc4
* d0>d1 : il faut revenir en arrière d'un pas
       suba.l  #2,a0
       suba.l  #2,a2
       add.l   d2,d7
       subq.w  #1,d6
* On passe maintenant à la couleur suivante
rc4:
       move.w  d6,(a1)          remise à jour de ColorTab
       adda.l  #4,a1
       subq.l  #1,d5 
       bra.s   rc2
rcEX:

* Inscription du nombre minimum d'iterations dans NMin(a5). C'est le 
* premier élément non nul de Histo
       move.l  Histo(pc),a0
       moveq   #-1,d0
min1   addq.w  #1,d0
       tst.w   (a0)+
       beq.s   min1
       move.w  d0,NMin(a5)  


Colorie:
* recoloriage, si d5 non nul,sinon on se contente de reformer Histo.
* Si on vient des lignes précédentes, d5 ne peut pas être nul,
* sinon il y aurait eu une division par 0 en 'rc2'

       move.l  Histo(pc),a4

* on se positionne en début de "Iter", aprés l'entête
       move.l  IterHandle(pc),d1
       moveq   #8,d2        en position 8
       move.l  #-1,d3       à partir du début
       move.l  dosbase(pc),a6
       jsr     Seek(a6)

       move.l  Tradu(pc),a3
       moveq   #0,d7
       move.l  dosbase(pc),a6
rcligne:
       move.l  IterHandle(pc),d1
       move.l  Ligne(pc),d2
       move.l  d2,a2
       move.w  NX(a5),d3
       ext.l   d3
       add.w   d3,d3
       jsr     Read(a6)
       moveq   #0,d6
rcpix:
       move.w  (a2)+,d1
       add.w   d1,d1
       move.w  0(a3,d1.w),d0
       tst.w   d5
       bne.s   rcp1
       add.w   #1,0(a4,d1.w)       reconstruction de Histo
       bra.s   rcp2
rcp1:  bsr     WritePix       
rcp2:  addq.w  #1,d6
       cmp.w   NX(a5),d6
       bne.s   rcpix
       addq.w  #1,d7
       cmp.w   NY(a5),d7
       bne.s   rcligne     
       
* Si d5=0 on revient en rc0
       tst.w   d5
       beq     rc0


fin:   
*       move.l  dosbase(pc),a6
*       jsr     IoErr(a6)

       move.l  intuibase(pc),a6
       move.l  Window(a5),a0
       jsr     -60(a6)             RefreshPointer

       move.l  ExecBase,a6
       move.l  MousePointer(pc),d0
       beq.s   fin0
       move.l  d0,a1
       move.l  #108,d0
       jsr     FreeMem(a6)

fin0:  move.l  Ligne(pc),d0
       beq.s   fin1
       move.l  d0,a1
       move.w  NX(a5),d0
       ext.l   d0
       asl     #1,d0
       cmp.b   #2,Iter(a5)
       bne.s   fina 
       asl     #5,d0
fina:  jsr     FreeMem(a6)

fin1:
       move.l  Histo(pc),d0
       beq.s   fin1a
       move.l  d0,a1
       move.w  NMax(a5),d0
       ext.l   d0
       addq.w  #1,d0
       add.l   d0,d0
       jsr     FreeMem(a6)
fin1a:
       move.l  PixTable(pc),d0
       beq.s   fin1b
       move.l  d0,a1
       move.w  NY(a5),d0
       ext.l   d0
       asl.l   #2,d0
       jsr     FreeMem(a6)
fin1b:
       move.l  Tradu(pc),d0
       beq.s   fin2
       move.l  d0,a1
       move.w  NMax(a5),d0
       ext.l   d0
       addq.w  #1,d0
       add.l   d0,d0
       jsr     FreeMem(a6)

fin2:
       lea     IterHandle(pc),a2
       move.l  (a2),d1
       beq.s   fin3
       move.l  dosbase(pc),a6
       jsr     Close(a6)
fin3:               
       move.l  #0,(a2)+
       move.l  #0,(a2)+
       move.l  #0,(a2)+
       move.l  #0,(a2)

   rts


* xxxxxxxxxxxxxxxxxxxxx Ouverture fichier, buffers... xxxxxxxxxxxxxx

OpenAll:



* Chargement du nouveau pointeur
       move.l ExecBase,a6
       move.l #65538,d1
       move.l #108,d0
       jsr    AllocMem(a6)
       PUSH   l,d0,MousePointer
       beq    OpenExit
       move.l d0,a0
       lea    pointer(pc),a1
       move.l #53,d1
ptr:   move.w (a1)+,(a0)+
       dbf    d1,ptr
       move.l Window(a5),a0 
       move.l d0,a1            MousePointer
       moveq  #25,d0
       moveq  #16,d1
       moveq  #0,d2
       moveq  #0,d3
       move.l intuibase(pc),a6
       jsr    -270(a6)         SetPointer         


* Ouverture d'un buffer 'Ligne' pour recevoir les résultats de calcul
* sur 1 ligne (2*NX octets)  si (Iter) différent de 2, sinon
* sur 32 lignes de NX mots, soit 64*NX octets. Le nombre de lignes
* est stocké dans NbLinBuffer
 
       lea NbLinBuffer(pc),a0
       move.w #1,(a0)
       cmp.b  #2,Iter(a5)
       bne.s  opa1
       move.w #32,(a0)
opa1:
       move.l ExecBase,a6
       moveq  #0,d0
       move.w NX(a5),d0
       asl    #1,d0
       cmp.b  #2,Iter(a5)
       bne.s  opa2
       asl    #5,d0
opa2:  move.l #$10001,d1
       jsr    AllocMem(a6)
       PUSH   l,d0,Ligne
       beq    OpenExit

* Ouverture d'une table de NY pointeurs (longs) pour Write/ReadPix
       moveq  #0,d0
       move.w NY(a5),d0
       asl.l  #2,d0
       move.l #$10001,d1
       jsr    AllocMem(a6)
       PUSH   l,d0,PixTable
       beq    OpenExit


* Ouverture d'un buffer 'Histo' pour l'histogramme
       moveq   #0,d0
       move.w  NMax(a5),d0
       addq.l  #1,d0
       add.l   d0,d0
       move.l  #$10001,d1
       movem.l d0/d1,-(sp)
       jsr     AllocMem(a6)
       PUSH    l,d0,Histo
       movem.l (sp)+,d0/d1
       beq.s   OpenExit

* Ouverture d'un buffer 'Tradu' pour la conversion des couleurs
       jsr    AllocMem(a6)
       PUSH   l,d0,Tradu
       beq.s  OpenExit
       

* Remplissage de 'Tradu'

       move.l ColTab(a5),a0
       addq.l #2,a0                a0 pointe vers la 1ère couleur
       move.l Tradu(pc),a1
       moveq  #0,d0
tra0:  cmp.w  2(a0),d0
       bls.s  tra1
       adda.l #4,a0
tra1:  move.w (a0),(a1)+
       addq.w #1,d0
       cmp.w  NMax(a5),d0
       bne.s  tra0
       move.w #0,(a1)


* Si Args(0)=(a5)='M' ou 'J', ouverture d'un fichier "iter" pour
* stocker les résultats de calcul, sauf si Iter(a5) est nul.
* Si (a5)='r' ou 'R', ouverture en lecture; dans ce cas, si le fichier
* n'existe pas, sortie d'échec (-> beq failure).

       move.l #1005,d2            OLDFILE
       move.b (a5),d7
       and.b  #$DF,d7
       cmp.b  #'R',d7
       beq.s  iter0
       move.l #1006,d2            NEWFILE
       move.b Iter(a5),d0
       bne.s  iter0
       moveq  #1,d0
       rts

iter0: moveq  #0,d1
       lea    IterHandle(pc),a2
       move.l d1,(a2)             efface le handle avant réouverture
       move.l dosbase(pc),a6
       lea    IterName(pc),a1
       move.l a1,d1
       cmp.b  #1,d0               Si (Iter)=1   -> "Iter" en RAM:
       beq.s  iter1
       move.l Work(a5),d1         sinonà l'adresse dans Args(0)+Work

iter1: jsr    Open(a6)
       move.l d0,(a2)

OpenExit:  rts


* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

* Routines d'écriture et de lecture d'un pixel à l'écran, inspirée de
* FastPix, de Scott B. Steiman (AC 4.11 p.23, 1989)
* On doit d'abord initialiser deux tables :
* -  PixTable contient les adresses du 1er octet de chacune des lignes du
*    dernier plan-bit (NY longs),
* -  EcartBp contient les écarts des adresses entre 2 plans-bits
*    consécutifs (depth-1 longs)
* L'appel à WritePix ou ReadPix se fera avec X(raster) et Y(raster) dans
* d6 et d7. La couleur sera transmise ou fournie dans d0 
* Attention : d1-d3 modifiés

InitPix:
       move.l   RastPort(pc),a0
       move.l   4(a0),a0           BitMap
       moveq    #0,d0
       move.b   5(a0),d0           depth
       subq.b   #1,d0
       lea      Depth1(pc),a1
       move.w   d0,(a1)+           Depth-1
       move.w   (a0),d2 
       ext.l    d2                 d2=BytesPerRow (long)

* a1 pointe vers le 1er EcartBp. L'adresse du dernier plan-bit est dans 
* a0+8+4*depth1 .  On remplit EcartBp.
       adda.l   #8,a0
       move.l   d0,d1
       asl.w    #2,d1
       adda.l   d1,a0
       move.l   (a0),a2            sauve l'adresse du dernier plan-bit
       subq.w   #1,d0
ip1:   move.l   (a0),d1
       sub.l    -(a0),d1           adresse plan (n) - adresse plan (n-1)
       move.l   d1,(a1)+
       dbf      d0,ip1

* Remplissage de PixTable
       move.l   PixTable(pc),a0
       move.l   a2,(a0)+
       move.w   NY(a5),d0
       subq.w   #2,d0
ip2:   add.l    d2,a2
       move.l   a2,(a0)+
       dbf      d0,ip2

       rts


Depth1:   dc.w   0         variables internes aux routines de pixels
EcartBp:  dcb.l  5,0

* - - - - - - - - - - - - - - - - - - - - -
* L'adresse de l'octet à modifier par WritePix dans le dernier plan-bit
* est PixTable+4*d7+d6/8. Les adresses dans les plans bits précédents
* s'obtiennent en retranchant les EcartBp successifs

WritePix:
       move.l    PixTable(pc),a1
       move.l    d7,d2
       asl.l     #2,d2
       move.l    0(a1,d2.w),a1     => 1er octet à modifier
       move.l    d6,d2
       asr.l     #3,d2             rang de l'octet dans la ligne
       adda.l    d2,a1

       move.l    d6,d1
       moveq     #7,d3
       and.w     d3,d1             masque
       sub.w     d1,d3
       move.w    Depth1(pc),d1
       lea       EcartBp(pc),a0
wp1:
       btst      d1,d0
       bne.s     wp2
       bclr      d3,(a1)           clears pixel
       sub.l     (a0)+,a1
       dbf       d1,wp1
       rts 
wp2:
       bset      d3,(a1)
       sub.l     (a0)+,a1
       dbf       d1,wp1
       rts

* ----------------------------------------------------------      
* L'adresse de l'octet à lire par ReadPix dans le dernier plan-bit
* est PixTable+4*d7+d6/8. Les adresses dans les plans bits précédents
* s'obtiennent en retranchant les EcartBp successifs

ReadPix:
       moveq     #0,d0
       move.l    PixTable(pc),a1
       move.l    d7,d2
       asl.l     #2,d2
       move.l    0(a1,d2.w),a1     => 1er octet à lire
       move.l    d6,d2
       asr.l     #3,d2             rang de l'octet dans la ligne
       adda.l    d2,a1

       move.l    d6,d1
       moveq     #7,d3
       and.w     d3,d1             masque
       sub.w     d1,d3
       move.w    Depth1(pc),d1
       lea       EcartBp(pc),a0
rp1:
       btst      d3,(a1)
       bne.s     rp2
       bclr      d1,d0            clears pixel
       sub.l     (a0)+,a1
       dbf       d1,rp1
       rts 
rp2:
       bset      d1,d0
       sub.l     (a0)+,a1
       dbf       d1,rp1
       rts


* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

*                |-----------------------------------|
*                |   CALCUL DU NOMBRE D'ITERATIONS   |
*                |-----------------------------------|
*
* X et Y chargés sur 64 bits en (a0,a1) et (a2,a3)
* On doit placer cx/2 et cy/2 dans (CX0,CX1) et (CY0,CY1)
* Sortie du nombre de tours en d0
* Tous les registres d, a0-3, a6 utilisés

Calcul:

     lea    Compteur(pc),a6
     move.w NMax(a5),(a6)
     lea    CX0(pc),a0

* Mandelbrot ou Julia ?
     cmp.b  #'J',(a5)
     bne.s  Mandel
     move.l XS(a5),d0           Julia
     move.l XS+4(a5),d1
     asr.l  #1,d0
     roxr.l #1,d1   
     move.l d0,(a0)+
     move.l d1,(a0)+  
     move.l YS(a5),d0
     move.l YS+4(a5),d1
     asr.l  #1,d0
     roxr.l #1,d1   
     move.l d0,(a0)+
     move.l d1,(a0)
     movea.l Xcur0(pc),a0  
     movea.l Xcur1(pc),a1  
     movea.l Ycur0(pc),a2  
     movea.l Ycur1(pc),a3  
     bra    BoucleCalcul

Mandel:
     move.l Xcur0(pc),d0    chargement de CX0,CX1 et CY0,CY1
     move.l Xcur1(pc),d1    et initialisations pour Mandelbrot
     move.l d1,a1
     asr.l  #1,d0
     roxr.l #1,d1
     move.l d0,(a0)+
     move.l d1,(a0)+
     move.l Ycur0(pc),d0
     move.l Ycur1(pc),d1
     move.l d0,a2
     move.l d1,a3
     asr.l  #1,d0
     roxr.l #1,d1
     move.l d0,(a0)+
     move.l d1,(a0)
     movea.l Xcur0(pc),a0

     cmp.b  #'M',(a5)
     beq    BoucleCalcul
     cmp.b  #'1',(a5)
     bne.s  Mandel2
     movea.l XS(a5),a0     Initialisation MandelIncomplet-1
     movea.l YS(a5),a2
     movea.l #0,a1
     movea.l #0,a3

Mandel2
     cmp.b   #'2',(a5)
     bne.s   Mandel3
     lea     CX0(pc),a0     MandelIncomplet-2
     move.l  YS(a5),d0
     move.l  YS+4(a5),d1
     asr.l   #1,d0
     roxr.l  #1,d1
     move.l  d0,(a0)+
     move.l  d1,(a0)+
     move.l  Ycur0(pc),d0
     move.l  Ycur1(pc),d1
     asr.l   #1,d0
     roxr.l  #1,d1
     move.l  d0,(a0)+
     move.l  d1,(a0)
     move.l  XS(a5),a0
     move.l  XS+4(a5),a1
     move.l  Xcur0(pc),a2
     move.l  Xcur1(pc),a3
     bra   BoucleCalcul

Mandel3
     cmp.b   #'3',(a5)
     bne.s   Mandel4
     lea     CX0(pc),a0        MandelIncomplet-3
     move.l  Xcur0(pc),d0
     move.l  Xcur1(pc),d1
     asr.l   #1,d0
     roxr.l  #1,d1
     move.l  d0,(a0)+
     move.l  d1,(a0)+
     move.l  YS(a5),d0
     move.l  YS+4(a5),d1
     asr.l   #1,d0
     roxr.l  #1,d1
     move.l  d0,(a0)+
     move.l  d1,(a0)
     move.l  XS(a5),a0
     move.l  XS+4(a5),a1
     move.l  Ycur0(pc),a2
     move.l  Ycur1(pc),a3
     bra.s   BoucleCalcul

Mandel4
     cmp.b   #'4',(a5)
     bne.s   Mandel5
     lea     CX0(pc),a0        MandelIncomplet-4
     move.l  XS(a5),d0
     move.l  XS+4(a5),d1
     asr.l   #1,d0
     roxr.l  #1,d1
     move.l  d0,(a0)+
     move.l  d1,(a0)+
     move.l  Ycur0(pc),d0
     move.l  Ycur1(pc),d1
     asr.l   #1,d0
     roxr.l  #1,d1
     move.l  d0,(a0)+
     move.l  d1,(a0)
     move.l  Xcur0(pc),a0
     move.l  Xcur1(pc),a1
     move.l  YS(a5),a2
     move.l  YS+4(a5),a3
     bra.s   BoucleCalcul

Mandel5
     cmp.b   #'5',(a5)
     bne.s   BoucleCalcul 
     lea     CX0(pc),a0           MandelIncomplet-5
     move.l  Xcur0(pc),d0
     move.l  Xcur1(pc),d1
     asr.l   #1,d0
     roxr.l  #1,d1
     move.l  d0,(a0)+
     move.l  d1,(a0)+
     move.l  XS(a5),d0
     move.l  XS+4(a5),d1
     asr.l   #1,d0
     roxr.l  #1,d1
     move.l  d0,(a0)+
     move.l  d1,(a0)
     move.l  YS(a5),a2
     move.l  YS+4(a5),a3
     move.l  Ycur0(pc),a0
     move.l  Ycur1(pc),a1

BoucleCalcul:
     cmpi.b #48,BitNb(a5)
     beq    BoucleCalcul48
     bls    BoucleCalcul32

BoucleCalcul64:

*======================  U = X**2 ============================ 
* Carré 64 bits de (a0,a1) sur (Destination)  (d4,d5)
* A partir de la décomposition en mots de 16 bits (a0,a1)=(N0,N1,N2,N3)
* le carré s'écrit      N0*N0 sur bits  1 à 32 (d4)
*                     2*N0*N1 sur bits 17 à 48 (d4.w et d5.haut)
*               N1*N1+2*N0*N2 sur bits 33 à 64 (d5)
*             2*N2*N1+2*N0*N3 sur bits 49 à 64 (d5.w)

XX:
      move.l a1,d5      N3 en place
      move.l a0,d1      N1 en place
      bpl.s  XOK
       neg.l  d5
       negx.l d1
XOK:  move.l d5,d2
      swap   d2         N2 en place
      move.l d1,d0
      swap   d0         N0 en place
      move.l #$FFFF,d4

* 2*N2*N1+2*N0*N3 sur bits 49 à 64 (d5.w)
      mulu   d0,d5
      swap   d5 
      and.l  d4,d5      N0*N3 dans d5.w
      add.l  d5,d5     2N0*N3
      move.w d2,d3
      mulu   d1,d3      N1*N2
      swap   d3
      and.l  d4,d3      N1*N2 dans d3.w
      add.l  d3,d5
      add.l  d3,d5

* N0*N0 sur bits  1 à 32 (d4)
      move.w d0,d4
      mulu   d4,d4       N0*N0 en place
     
      mulu   d0,d2       N0*N2 dans d2
      mulu   d1,d0       N1*N0 dans d0
      mulu   d1,d1       N1*N1 dans d1

* N1*N1+2*N0*N2 sur bits 33 à 64 (d5)
      moveq  #0,d3
      add.l  d2,d5        en place
      addx.l d3,d4        avec la retenue
      add.l  d2,d5       et le 2ème
      addx.l d3,d4        avec la retenue
      add.l  d1,d5        N1*N1 en place
      addx.l d3,d4         avec la retenue

* d1 et d2 sont maintenant disponibles
      move.l #$FFFF,d1
      
* 2*N0*N1 sur bits 17 à 48 (d4.w et d5.haut)
      add.l  d0,d0         2*N0*N1
      bcc.s  sq0
      add.l  #$10000,d4    retenue
sq0:  move.l d0,d2
      and.l  d1,d2
      swap   d2            bits bas de 2*N0*N1 dans d2.haut
      add.l  d2,d5
      addx.l d3,d4   
      swap   d0
      and.l  d1,d0         bits hauts de 2*N0*N1 dans d0.w
      add.l  d0,d4         pas de problème de retenue à craindre

*======================  V = Y**2 ============================ 
* Carré 64 bits de (a2,a3) sur (Destination)  (d6,d7)
* A partir de la décomposition en mots de 16 bits (a0,a1)=(N0,N1,N2,N3)
* le carré s'écrit      N0*N0 sur bits  1 à 32 (d6)
*                     2*N0*N1 sur bits 17 à 48 (d6.w et d7.haut)
*               N1*N1+2*N0*N2 sur bits 33 à 64 (d7)
*             2*N2*N1+2*N0*N3 sur bits 49 à 64 (d7.w)

YY:
      move.l a3,d7      N3 en place
      move.l a2,d1      N1 en place
      bpl.s  YOK
       neg.l  d7
       negx.l d1
YOK:  move.l d7,d2
      swap   d2         N2 en place
      move.l d1,d0
      swap   d0         N0 en place
      move.l #$FFFF,d6

* 2*N2*N1+2*N0*N3 sur bits 49 à 64 (d7.w)
      mulu   d0,d7
      swap   d7 
      and.l  d6,d7      N0*N3 dans d7.w
      add.l  d7,d7     2N0*N3
      move.w d2,d3
      mulu   d1,d3      N1*N2
      swap   d3
      and.l  d6,d3      N1*N2 dans d3.w
      add.l  d3,d7
      add.l  d3,d7

* N0*N0 sur bits  1 à 32 (d6)
      move.w d0,d6
      mulu   d6,d6       N0*N0 en place
     
      mulu   d0,d2       N0*N2 dans d2
      mulu   d1,d0       N1*N0 dans d0
      mulu   d1,d1       N1*N1 dans d1

* N1*N1+2*N0*N2 sur bits 33 à 64 (d7)
      moveq  #0,d3
      add.l  d2,d7        en place
      addx.l d3,d6        avec la retenue
      add.l  d2,d7       et le 2ème
      addx.l d3,d6        avec la retenue
      add.l  d1,d7        N1*N1 en place
      addx.l d3,d6         avec la retenue

* d1 et d2 sont maintenant disponibles
      move.l #$FFFF,d1
      
* 2*N0*N1 sur bits 17 à 48 (d6.w et d7.haut)
      add.l  d0,d0         2*N0*N1
      bcc.s  sq1
      add.l  #$10000,d6    retenue
sq1:  move.l d0,d2
      and.l  d1,d2
      swap   d2            bits bas de 2*N0*N1 dans d2.haut
      add.l  d2,d7
      addx.l d3,d6   
      swap   d0
      and.l  d1,d0         bits hauts de 2*N0*N1 dans d0.w
      add.l  d0,d6         pas de problème de retenue à craindre

*================ S = U+V dans (d6,d7) =========================
*                     U-V dans (d4,d5)
* Sortie si S trop grand

UPV:
      move.l d6,d0
      move.l d7,d1
      add.l  d5,d7
      addx.l d4,d6
      sub.l  d1,d5
      subx.l d0,d4

      cmpi.l #$40000000,d6
      bcc    iterexit
      sub.w   #1,(a6)
      beq    iterexit


*================== |X+Y| dans (d0,d1) ===========================

XPY:
      move.l  a0,d0
      move.l  a1,d1
      move.l  a2,d2
      add.l   a3,d1
      addx.l  d2,d0
      bge.s   it1
      neg.l   d1
      negx.l  d0
it1:

*=================== Mise à jour de X =============================
*                    X = U-V + CX
* Travail dans (d4,d5) (U-V au départ) puis transfert en (a0,a1)

      add.l   d5,d5
      addx.l  d4,d4
      move.l  CX0(pc),d2       CX/2, sur 64 bits
      add.l   CX1(pc),d5   
      addx.l  d2,d4
      bvs     iterexit
* Repositionnement du point décimal
      add.l   d5,d5
      addx.l  d4,d4
      bvs     iterexit
* Transfert
      move.l  d5,a1
      move.l  d4,a0
      
*==================== |X+Y|**2 dans (d4,d5) ==========================
*                   ( (d0,d1) au départ)
* Carré 64 bits de (d0,d1) sur (Destination)  (d4,d5)
* A partir de la décomposition en mots de 16 bits (d0,d1)=(N0,N1,N2,N3)
* le carré s'écrit      N0*N0 sur bits  1 à 32 (d4)
*                     2*N0*N1 sur bits 17 à 48 (d4.w et d5.haut)
*               N1*N1+2*N0*N2 sur bits 33 à 64 (d5)
*             2*N2*N1+2*N0*N3 sur bits 49 à 64 (d5.w)

XY2:
      move.l d1,d5      N3 en place
      move.l d5,d2
      swap   d2         N2 en place
      move.l d0,d1      N1 en place
      move.l d1,d0
      swap   d0         N0 en place
      move.l #$FFFF,d4

* 2*N2*N1+2*N0*N3 sur bits 49 à 64 (d5.w)
      mulu   d0,d5
      swap   d5 
      and.l  d4,d5      N0*N3 dans d5.w
      add.l  d5,d5     2N0*N3
      move.w d2,d3
      mulu   d1,d3      N1*N2
      swap   d3
      and.l  d4,d3      N1*N2 dans d3.w
      add.l  d3,d5
      add.l  d3,d5

* N0*N0 sur bits  1 à 32 (d4)
      move.w d0,d4
      mulu   d4,d4       N0*N0 en place
     
      mulu   d0,d2       N0*N2 dans d2
      mulu   d1,d0       N1*N0 dans d0
      mulu   d1,d1       N1*N1 dans d1

* N1*N1+2*N0*N2 sur bits 33 à 64 (d5)
      moveq  #0,d3
      add.l  d2,d5        en place
      addx.l d3,d4        avec la retenue
      add.l  d2,d5       et le 2ème
      addx.l d3,d4        avec la retenue
      add.l  d1,d5        N1*N1 en place
      addx.l d3,d4         avec la retenue

* d1 et d2 sont maintenant disponibles
      move.l #$FFFF,d1
      
* 2*N0*N1 sur bits 17 à 48 (d4.w et d5.haut)
      add.l  d0,d0         2*N0*N1
      bcc.s  sq2
      add.l  #$10000,d4    retenue
sq2:  move.l d0,d2
      and.l  d1,d2
      swap   d2            bits bas de 2*N0*N1 dans d2.haut
      add.l  d2,d5
      addx.l d3,d4   
      swap   d0
      and.l  d1,d0         bits hauts de 2*N0*N1 dans d0.w
      add.l  d0,d4         pas de problème de retenue à craindre

*=================== Mise à jour de Y ============================
*                   Y = |X+Y|**2 -S + CY
* Travail dans (d4,d5)

NEWY:
      sub.l  d7,d5
      subx.l d6,d4         pas de probleme de débordement

      add.l   d5,d5
      addx.l  d4,d4
      move.l  CY0(pc),d2       CY/2, sur 64 bits
      add.l   CY1(pc),d5
      addx.l  d2,d4
      bvs.s     iterexit
* Repositionnement du point décimal
      add.l   d5,d5
      addx.l  d4,d4
      bvs.s   iterexit
* Transfert
      move.l  d5,a3
      move.l  d4,a2

*======================= Un tour de plus ====================

      bra     BoucleCalcul64

*============================================================

* C'est fini!
iterexit:
      move.w  NMax(a5),d0
      sub.w   (a6),d0
      rts



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

BoucleCalcul48:

*======================  U = X**2 ============================ 
* Carré 48 bits de (a0,a1) sur (Destination)  (d4,d5)
* A partir de la décomposition en mots de 16 bits (a0,a1)=(N0,N1,N2,N3)
* le carré s'écrit      N0*N0 sur bits  1 à 32 (d4)
*                     2*N0*N1 sur bits 17 à 48 (d4.w et d5.haut)
*               N1*N1+2*N0*N2 sur bits 33 à 64 (d5)

      move.l a1,d2      (N3 en place)
      move.l a0,d1      N1 en place
      bpl.s  XOK1
       neg.l  d2
       negx.l d1
XOK1: 
      swap   d2         N2 en place
      move.l d1,d0
      swap   d0         N0 en place
      move.l #$FFFF,d4


* N0*N0 sur bits  1 à 32 (d4)
      move.w d0,d4
      mulu   d4,d4       N0*N0 en place
     
      mulu   d0,d2       N0*N2 dans d2
      mulu   d1,d0       N1*N0 dans d0
      mulu   d1,d1       N1*N1 dans d1

* N1*N1+2*N0*N2 sur bits 33 à 64 (d5)
      moveq  #0,d3
      move.l d2,d5       N0*N2 en place (une fois)
      add.l  d2,d5           et 2 fois
      addx.l d3,d4        avec la retenue
      add.l  d1,d5        N1*N1 en place
      addx.l d3,d4         avec la retenue

* d1 et d2 sont maintenant disponibles
      move.l #$FFFF,d1
      
* 2*N0*N1 sur bits 17 à 48 (d4.w et d5.haut)
      add.l  d0,d0         2*N0*N1
      bcc.s  sqx1
      add.l  #$10000,d4    retenue
sqx1:  move.l d0,d2
      and.l  d1,d2
      swap   d2            bits bas de 2*N0*N1 dans d2.haut
      add.l  d2,d5
      addx.l d3,d4   
      swap   d0
      and.l  d1,d0         bits hauts de 2*N0*N1 dans d0.w
      add.l  d0,d4         pas de problème de retenue à craindre

*======================  V = Y**2 ============================ 
* Carré 48 bits de (a2,a3) sur (Destination)  (d6,d7)
* A partir de la décomposition en mots de 16 bits (a0,a1)=(N0,N1,N2,N3)
* le carré s'écrit      N0*N0 sur bits  1 à 32 (d6)
*                     2*N0*N1 sur bits 17 à 48 (d6.w et d7.haut)
*               N1*N1+2*N0*N2 sur bits 33 à 64 (d7)


      move.l a3,d2  
      move.l a2,d1      N1 en place
      bpl.s  YOK1
       neg.l  d2
       negx.l d1
YOK1: 
      swap   d2         N2 en place
      move.l d1,d0
      swap   d0         N0 en place
      move.l #$FFFF,d6


* N0*N0 sur bits  1 à 32 (d6)
      move.w d0,d6
      mulu   d6,d6       N0*N0 en place
     
      mulu   d0,d2       N0*N2 dans d2
      mulu   d1,d0       N1*N0 dans d0
      mulu   d1,d1       N1*N1 dans d1

* N1*N1+2*N0*N2 sur bits 33 à 64 (d7)
      moveq  #0,d3
      move.l d2,d7        en place
      add.l  d2,d7       et le 2ème
      addx.l d3,d6        avec la retenue
      add.l  d1,d7        N1*N1 en place
      addx.l d3,d6         avec la retenue

* d1 et d2 sont maintenant disponibles
      move.l #$FFFF,d1
      
* 2*N0*N1 sur bits 17 à 48 (d6.w et d7.haut)
      add.l  d0,d0         2*N0*N1
      bcc.s  sqy1
      add.l  #$10000,d6    retenue
sqy1:  move.l d0,d2
      and.l  d1,d2
      swap   d2            bits bas de 2*N0*N1 dans d2.haut
      add.l  d2,d7
      addx.l d3,d6   
      swap   d0
      and.l  d1,d0         bits hauts de 2*N0*N1 dans d0.w
      add.l  d0,d6         pas de problème de retenue à craindre

*================ S = U+V dans (d6,d7) =========================
*                     U-V dans (d4,d5)
* Sortie si S trop grand

      move.l d6,d0
      move.l d7,d1
      add.l  d5,d7
      addx.l d4,d6
      sub.l  d1,d5
      subx.l d0,d4

      cmpi.l #$40000000,d6
      bcc    iterexit48
      sub.w   #1,(a6)
      beq    iterexit48


*================== |X+Y| dans (d0,d1) ===========================


      move.l  a0,d0
      move.l  a1,d1
      move.l  a2,d2
      add.l   a3,d1
      addx.l  d2,d0
      bge.s   it11
      neg.l   d1
      negx.l  d0
it11:

*=================== Mise à jour de X =============================
*                    X = U-V + CX
* Travail dans (d4,d5) (U-V au départ) puis transfert en (a0,a1)

      add.l   d5,d5
      addx.l  d4,d4
      move.l  CX0(pc),d2       CX/2, sur 64 bits
      add.l   CX1(pc),d5   
      addx.l  d2,d4
      bvs.s     iterexit48
* Repositionnement du point décimal
      add.l   d5,d5
      addx.l  d4,d4
      bvs.s     iterexit48
* Transfert
      move.l  d5,a1
      move.l  d4,a0
      
*==================== |X+Y|**2 dans (d4,d5) ==========================
*                   ( (d0,d1) au départ)
* Carré 48 bits de (d0,d1) sur (Destination)  (d4,d5)
* A partir de la décomposition en mots de 16 bits (d0,d1)=(N0,N1,N2,N3)
* le carré s'écrit      N0*N0 sur bits  1 à 32 (d4)
*                     2*N0*N1 sur bits 17 à 48 (d4.w et d5.haut)
*               N1*N1+2*N0*N2 sur bits 33 à 64 (d5)

      move.l d1,d2
      swap   d2         N2 en place
      move.l d0,d1      N1 en place
      swap   d0         N0 en place


* N0*N0 sur bits  1 à 32 (d4)
      move.w d0,d4
      mulu   d4,d4       N0*N0 en place
     
      mulu   d0,d2       N0*N2 dans d2
      mulu   d1,d0       N1*N0 dans d0
      mulu   d1,d1       N1*N1 dans d1

* N1*N1+2*N0*N2 sur bits 33 à 64 (d5)
      moveq  #0,d3
      move.l d2,d5        en place
      add.l  d2,d5       et le 2ème
      addx.l d3,d4        avec la retenue
      add.l  d1,d5        N1*N1 en place
      addx.l d3,d4         avec la retenue

* d1 et d2 sont maintenant disponibles
      move.l #$FFFF,d1
      
* 2*N0*N1 sur bits 17 à 48 (d4.w et d5.haut)
      add.l  d0,d0         2*N0*N1
      bcc.s  sqz2
      add.l  #$10000,d4    retenue
sqz2:  move.l d0,d2
      and.l  d1,d2
      swap   d2            bits bas de 2*N0*N1 dans d2.haut
      add.l  d2,d5
      addx.l d3,d4   
      swap   d0
      and.l  d1,d0         bits hauts de 2*N0*N1 dans d0.w
      add.l  d0,d4         pas de problème de retenue à craindre

*=================== Mise à jour de Y ============================
*                   Y = |X+Y|**2 -S + CY
* Travail dans (d4,d5)

      sub.l  d7,d5
      subx.l d6,d4         pas de probleme de débordement

      add.l   d5,d5
      addx.l  d4,d4
      move.l  CY0(pc),d2       CY/2, sur 64 bits
      add.l   CY1(pc),d5
      addx.l  d2,d4
      bvs.s     iterexit48
* Repositionnement du point décimal
      add.l   d5,d5
      addx.l  d4,d4
      bvs.s   iterexit48
* Transfert
      move.l  d5,a3
      move.l  d4,a2

*======================= Un tour de plus ====================

      bra     BoucleCalcul48

*============================================================

* C'est fini!
iterexit48:
      move.w  NMax(a5),d0
      sub.w   (a6),d0
      rts

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


BoucleCalcul32:

* On va utiliser la recette de MandelVroom pour raccourcir les calculs dans
* l'ensemble principal. Selon qu'on sort avec N=NMax ou non, on met
* le flag SetFlag(a5) à 1 ou 0. Si ce flag est à 1 (c.à.d. si on vient de
* rentrer dans l'ensemble), on mémorise les derniers x2+y2 dans une table
* (les 16 derniers dans MandelVroom) et quand la table est pleine, on 
* examine si le dernier nombre n'était pas déjà dans la table. Si oui, on
* est dans l'ensemble principal. 
* On se sert des 42 octets de la zone "TableTraces" pour le stockage

table=20
	move.w   SetFlag(a5),d7
        lea      TableTraces(pc),a1
        move.w   #table,d5


Boucle32:

*======================  U = X**2 ============================ 
* Carré 32 bits de (a0) sur (Destination)  (d4)
* A partir de la décomposition en mots de 16 bits (a0)=(N0,N1)
* le carré s'écrit      N0*N0 sur bits  1 à 32 (d4)
*                     2*N0*N1 sur bits 17 à 32 (d4.w)

      move.l a0,d1      N1 en place
      bpl.s  XOK2
      neg.l  d1
XOK2: 
      move.l d1,d0
      swap   d0         N0 en place
      move.l #$FFFF,d4


* N0*N0 sur bits  1 à 32 (d4)
      move.w d0,d4
      mulu   d4,d4       N0*N0 en place
     
      mulu   d1,d0       N1*N0 dans d0
* d1 et d2 sont maintenant disponibles
      move.l #$FFFF,d1
      
* 2*N0*N1 sur bits 17 à 32 (d4.w)
      add.l  d0,d0         2*N0*N1
      bcc.s  sqx2
      add.l  #$10000,d4    retenue
sqx2: 
      swap   d0
      and.l  d1,d0         bits hauts de 2*N0*N1 dans d0.w
      add.l  d0,d4         pas de problème de retenue à craindre

*======================  V = Y**2 ============================ 
* Carré 32 bits de (a2) sur (Destination)  (d6)
* A partir de la décomposition en mots de 16 bits (a2)=(N0,N1)
* le carré s'écrit      N0*N0 sur bits  1 à 32 (d6)
*                     2*N0*N1 sur bits 17 à 48 (d6.w)


      move.l a2,d1      N1 en place
      bpl.s  YOK2
       neg.l  d1
YOK2: 
      move.l d1,d0
      swap   d0         N0 en place
      move.l #$FFFF,d6


* N0*N0 sur bits  1 à 32 (d6)
      move.w d0,d6
      mulu   d6,d6       N0*N0 en place
     
      mulu   d1,d0       N1*N0 dans d0
* d1 et d2 sont maintenant disponibles
      move.l #$FFFF,d1
      
* 2*N0*N1 sur bits 17 à 32 (d6.w)
      add.l  d0,d0         2*N0*N1
      bcc.s  sqy2
      add.l  #$10000,d6    retenue
sqy2:  
      swap   d0
      and.l  d1,d0         bits hauts de 2*N0*N1 dans d0.w
      add.l  d0,d6         pas de problème de retenue à craindre

*================ S = U+V dans (d6) =========================
*                     U-V dans (d4)
* Sortie si S trop grand

      move.l d6,d0
      add.l d4,d6
      sub.l d0,d4

* remplissage de la table des traces (le raccourci MandelVroom)
      tst.w  d7         seulement si SetFlag (d7)=1 
      beq.s  debor32
*      bra.s  debor32    ** à réactiver pour calcul complet
      move.l d6,(a1)+   stokage d6 dans la table
      subq.w #1,d5
      bne.s  debor32    on examine la table au bout de 'table' tours
      move.w #table-2,d5
      suba.l #4,a1
examen:
      cmp.l  -(a1),d6
      beq.s  SetOK      on doit être dans le Set
      dbra   d5,examen      

      move.w #table,d5  on n'a pas trouvé. On recharge d5; a1 doit
      bra.s  debor32      pointer vers Histo à nouveau

SetOK: move.w NMax(a5),d0
       bra.s  ex31
 
* y-a-t-il eu débordement?
debor32:
      cmpi.l  #$40000000,d6
      bcc.s   iterexit32
      sub.w   #1,(a6)
      beq.s   iterexit32



*================== |X+Y| dans (d0) ===========================


      move.l  a0,d0
      add.l   a2,d0
      bge.s   it12
      neg.l   d0
it12:

*=================== Mise à jour de X =============================
*                    X = U-V + CX
* Travail dans (d4) (U-V au départ) puis transfert en (a0)

      add.l   d4,d4
      add.l   CX0(pc),d4
      bvs.s   iterexit32
* Repositionnement du point décimal
      add.l   d4,d4
      bvs.s   iterexit32
* Transfert
      move.l  d4,a0
      
*==================== |X+Y|**2 dans (d4,d5) ==========================
*                   ( (d0) au départ)
* Carré 32 bits de (d0) sur (Destination)  (d4)
* A partir de la décomposition en mots de 16 bits (d0)=(N0,N1)
* le carré s'écrit      N0*N0 sur bits  1 à 32 (d4)
*                     2*N0*N1 sur bits 17 à 48 (d4.w)

      move.l d0,d1      N1 en place
      swap   d0         N0 en place


* N0*N0 sur bits  1 à 32 (d4)
      move.w d0,d4
      mulu   d4,d4       N0*N0 en place
     
      mulu   d1,d0       N1*N0 dans d0
* d1 et d2 sont maintenant disponibles
      move.l #$FFFF,d1

      
* 2*N0*N1 sur bits 17 à 32 (d4.w)
      add.l  d0,d0         2*N0*N1
      bcc.s  sqz3
      add.l  #$10000,d4    retenue
sqz3: 
      swap   d0
      and.l  d1,d0         bits hauts de 2*N0*N1 dans d0.w
      add.l  d0,d4         pas de problème de retenue à craindre

*=================== Mise à jour de Y ============================
*                   Y = |X+Y|**2 -S + CY
* Travail dans (d4)

      sub.l  d6,d4  

      add.l   d4,d4
      add.l   CY0(pc),d4
      bvs.s   iterexit32
* Repositionnement du point décimal
      add.l   d4,d4
      bvs.s   iterexit32
* Transfert
      move.l  d4,a2

*======================= Un tour de plus ====================

      bra     Boucle32

*============================================================

* C'est fini!
iterexit32:
      move.w  NMax(a5),d0
      sub.w   (a6),d0
      move.w  #0,SetFlag(a5)
      cmp.w   NMax(a5),d0
      bne.s   ex32
ex31      move.w  #1,SetFlag(a5)
ex32      rts


***************** Division (64+1 bits) *************

* a0 pointe sur le dividende (sur 16+64 bits, mais le premier mot étant
*    0 ou 1).  Ce dividende est non signé
* a1 pointe sur le quotient
* le diviseur doit être dans d0

Div65:
      move.l (a0)+,d1
      moveq  #3,d2
divi: divu   d0,d1
      move.w d1,(a1)+        16 bits pour le quotient
      move.w (a0)+,d1        chargement des 16 bits suivants du dividende
      dbf    d2,divi
      rts

******************  (a3) = (a2) + (a1)*P(d1)/Q(d0) *****************
* Calculs sur 64 bits ; (a1) positif non signé
* P dans d1 (16 bits) et Q dans d0 (16 bits)
* Utilise un espace de travail  de 5+4 mots (a5+Work2)


MulDivAdd:
      bsr.s     MulDiv
      bsr.s     Add
      rts
MulDivSub:
      bsr.s   MulDiv
      bsr.s   Sub
      rts

MulDiv:
      move.l  a5,a0
      adda.l  #Work2,a0         pointe sur espace de travail
      moveq   #4,d2            qu'on va d'abord vider
mu1:  move.w  #0,(a0)+ 
      dbf     d2,mu1

* (a1)*P dans l'espace de travail
      suba.l  #2,a0            a0 => dernier mot du résultat
      adda.l  #8,a1            a1 => juste aprés le multiplicande
      moveq   #3,d2            pour 4 boucles
mu2:  move.w  -(a1),d3
      mulu    d1,d3
      suba.l  #2,a0
      add.l   d3,(a0)
      dbf     d2,mu2
* a0 pointe maintenant sur Work2. On va ranger P*Z1/Q dans Work2+10
      move.l  a0,a1  
      adda.l  #10,a1
      bsr.s   Div65
      rts

* a1 pointe juste aprés le quotient
Add:
      move.l  (a2)+,(a3)+
      move.l  (a2)+,(a3)+
      move.w  #0,ccr
      addx.l  -(a1),-(a3)  
      addx.l  -(a1),-(a3)  
      rts

Sub:
      move.l  (a2)+,(a3)+
      move.l  (a2)+,(a3)+
      move.w  #0,ccr
      subx.l  -(a1),-(a3)  
      subx.l  -(a1),-(a3)  
      rts


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

* Variables internes à l'iteration
CX0:   dc.l $7f0000              X courant, divisé par 2
CX1:   dc.l 0
CY0:   dc.l $0F0000              Y courant, divisé par 2
CY1:   dc.l 0
Xcur0: dc.l 0                    X courant
Xcur1: dc.l 0
Ycur0: dc.l 0                    Y courant
Ycur1: dc.l 0
IterHandle: dc.l 0                                    NE PAS
Ligne:      dc.l 0               => buffer "ligne"    TOUCHER
Histo:      dc.l 0                                    CES QUATRE
Tradu:      dc.l 0                                    LIGNES
PixTable:   dc.l 0               => table de pointeurs pour R/WPixel
MousePointer dc.l 0
NbLinBuffer  dc.w 1
Compteur:  dc.w 1000
TableTraces: dcb.l 21,0

X1=4        8
Y1=X1+8     8
X2=Y1+8
Y2=X2+8
NX=Y2+8
NY=NX+2
XR1=NY+2
YR1=XR1+2
XR2=YR1+2
DX=XR2+2    8
NX0=DX
NY0=NX0+2
DY=DX+8     8
XS=DY+8
YS=XS+8
Window=YS+8
NMax=Window+4
NbCol=NMax+2
ColTab=NbCol+2
Work=ColTab+4
CompteLignes=Work+4
SetFlag=CompteLignes+4
Work2=Work+5*2
X1J=Work2+9*2
Y1J=X1J+8
X2J=Y1J+8
Y2J=X2J+8
Total=Y2J+8
NMin=Work
BitNb=1        64 pour 64 bits, 48 pour 48 bits, 32 pour 32 bits
ModeCalcul=2    
Iter=3

intuitxt0:  dc.l  $1000,$40008,0,0,0
intuitxt1:  dc.l  $1000,$60003,0,0,0
intuitxt2:  dc.l  $1000,$60003,0,0,0
IterName:   dc.b  'ram:iter',0
text0:      dc.b  'Stopper le calcul?',0
text1:      dc.b  'OUI',0
text2:      dc.b  'NON',0
         even

* Sprite data pour le pointeur (27x2=54 mots=108 octets)
pointer: dc.w 0,0
         dc.w $8661,$8661
         dc.w $7c3e,$7c3e
         dc.w $0080,$0080
         dc.w $0080,$0080
         dc.w $0080,$0080
         dc.w $0080,$0080
         dc.w $0380,$0380
         dc.w 0,0
         dc.w $2004,$2004
         dc.w $1818,$1818
         dc.w $0ff0,$0ff0
         dc.w 0,0
         dc.w 0,0
         dc.w 0,$ee
         dc.w 0,$8a
         dc.w 0,$8e
         dc.w 0,$54ea
         dc.w 0,$40
         dc.w $804,0
         dc.w 0,4
         dc.w 0,$ab6E
         dc.w 0,$aa54
         dc.w 0,$ab54
         dc.w 0,$aa54
         dc.w 0,$6b56
         dc.w 0,0

         

fini:
  end