Fractale de Newton

' **********************************************************************
' Fractale de Newton
' **********************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image en pixels
dim PicHeight% : PicHeight% = 480

dim x0         : x0 = 0             : ' Coord. du centre de l'image
dim y0         : y0 = 0
dim Eps        : Eps = 0.000001     : ' Precision requise
dim MaxIter%   : MaxIter% = 100     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 1       : ' Facteur de zoom
dim ColorFact  : ColorFact = -2     : ' Facteur de coloration
dim CstVal     : CstVal = 0.9       : ' Luminosite HSV (constante)

' Variables supplementaires

dim HalfPicWidth%  : HalfPicWidth%  = PicWidth% / 2
dim HalfPicHeight% : HalfPicHeight% = PicHeight% / 2

dim R%, G%, B%  : ' Parametres de coloration
dim AbsColor    : ' abs(ColorFact)
dim ScaleFact   : ' Facteur d'echelle = distance entre 2 pixels
dim Nx%, Ny%    : ' Coordonnees d'un point (pixels)
dim xt, yt      : ' Coordonnees d'un point (algebriques)
dim f_x, f_y    : ' Fonction complexe f(z)
dim fp_x, fp_y  : ' Derivee f'(z)

' ----------------------------------------------------------------------
' Description des objets
' ----------------------------------------------------------------------

' Fenetre principale

left    0, 50
top     0, 50
width   0, PicWidth% + 70
height  0, PicHeight% + 120
caption 0, "Fractale de Newton ... Veuillez patienter."

' Image

picture 1
left    1, 30
top     1, 40
width   1, PicWidth%
height  1, PicHeight%

2d_target_is 1

' ----------------------------------------------------------------------
' Programme principal
' ----------------------------------------------------------------------

ColorFact = 0.01 * ColorFact
AbsColor  = abs(ColorFact)
ScaleFact = 4 / (PicHeight% * ZoomFact)

for Ny% = 0 to PicHeight% - 1
 yt = y0 - ScaleFact * (Ny% - HalfPicHeight%)
 for Nx% = 0 to PicWidth% - 1
   xt = x0 + ScaleFact * (Nx% - HalfPicWidth%)
   Newton(xt, yt)
   2d_pen_color R%, G%, B%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

file_save 1, "newton.bmp"

caption 0, "Fractale de Newton ... Terminé."

end

' ----------------------------------------------------------------------
' Sous-programmes
' ----------------------------------------------------------------------

sub Func(x, y)
' Calcul de la fonction complexe f(z) = (f_x, f_y)
' et de sa derivee f'(z) = (fp_x, fp_y)
' Ici f(z) = x^3 - 1 et f'(z) = 3 x^2

 dim_local x2, y2, x3, y3
 
 if x = 0 and y = 0
   
   f_x = 10
   f_y = 10
   
 else  
   
   x2 = x * x : x3 = x2 * x
   y2 = y * y : y3 = y2 * y
 
   f_x = x3 - 3 * x * y2 - 1
   f_y = 3 * x2 * y - y3
 
   fp_x = 3 * (x2 - y2)
   fp_y = 6 * x * y
 
 end_if    
end_sub

sub Newton(xt, yt)
' Iteration de la fonction complexe au point (xt, yt)

 dim_local Iter%, mf2, mq2, x, y, zn_x, zn_y, d, q_x, q_y

 zn_x = xt
 zn_y = yt
 
 Iter% = 0
 mf2 = 1
 mq2 = 1

 while Iter% < MaxIter% and mf2 > Eps and mq2 > Eps

   x = zn_x
   y = zn_y

   Func(x, y)

   d = fp_x * fp_x + fp_y * fp_y
 
   q_x = (f_x * fp_x + f_y * fp_y) / d  : ' q = f(z) / f'(z)  
   q_y = (f_y * fp_x - f_x * fp_y) / d
 
   zn_x = x - q_x
   zn_y = y - q_y

   mf2 = f_x * f_x + f_y * f_y
   mq2 = q_x * q_x + q_y * q_y
   
   Iter% = Iter% + 1

 end_while
 
 if Iter% > MaxIter%  : ' Pas de convergence
   R% = 0
   G% = 0
   B% = 0
 else
   MdbCol(Iter%)
 end_if
end_sub

sub MdbCol(Iter%)
' Determine la teinte (Hue, H) et la Saturation (S)
' d'apres le nombre d'iterations

 dim_local Q, Angle, Radius, H, S, V
 
 Q = log(Iter%) * AbsColor

 if Q < 0.5
   Q = 1 - 1.5 * Q
   Angle = 1 - Q
 else
   Q = 1.5 * Q - 0.5
   Angle = Q
 end_if

 Radius = sqr(Q)

' Si ColorFact > 0, assombrir une bande sur 2

 V = CstVal

 if (ColorFact > 0) and (odd(Iter%) > 0)
   V = 0.85 * V
   Radius = 0.667 * Radius
 end_if

 H = Angle * 10
 H = H - int(H)
 H = H * 360

 S = Radius - int(Radius)

' Convertir HSV en RGB

 HSVtoRGB(H, S, V)
end_sub

sub HSVtoRGB(H, S, V)
' Conversion HSV --> RGB. Resultats dans R%, G%, B%

 dim_local II%, ZZ, FF, PP, QQ, TT, RR, GG, BB
 
 if S = 0

   R% = int(V * 255)
   G% = R%
   B% = R%

 else

   ZZ  = H / 60
   II% = int(ZZ)
   FF  = ZZ - int(ZZ)
   PP  = V * (1 - S)
   QQ  = V * (1 - S * FF)
   TT  = V * (1 - S * (1 - FF))

   select II%
     case 0
       RR = V  : GG = TT : BB = PP
     case 1
       RR = QQ : GG = V  : BB = PP
     case 2
       RR = PP : GG = V  : BB = TT
     case 3
       RR = PP : GG = QQ : BB = V
     case 4
       RR = TT : GG = PP : BB = V
     case 5
       RR = V  : GG = PP : BB = QQ
   end_select

   R% = int(RR * 255)
   G% = int(GG * 255)
   B% = int(BB * 255)

 end_if
end_sub

' **************************************************************
' Fractale de Newton : Fonction z^p - 1 (p entier)
' **************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image
dim PicHeight% : PicHeight% = 480   : '   en pixels
dim p%         : p% = 5             : ' Exposant (z^p - 1)
dim x0         : x0 = 0             : ' Coordonnees du centre
dim y0         : y0 = 0             : '   de l'image
dim Eps        : Eps = 0.000001     : ' Precision requise
dim MaxIter%   : MaxIter% = 100     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 1       : ' Facteur de zoom
dim ColorFact  : ColorFact = -2     : ' Facteur de coloration
dim CstVal     : CstVal = 0.9       : ' Luminosite HSV constante

' Variables supplementaires

dim HalfPicWidth%  : HalfPicWidth%  = PicWidth% / 2
dim HalfPicHeight% : HalfPicHeight% = PicHeight% / 2
dim Eps2           : Eps2 = Eps * Eps

dim R%, G%, B%  : ' Parametres de coloration
dim AbsColor    : ' abs(ColorFact)
dim ScaleFact   : ' Facteur d'echelle = distance entre 2 pixels
dim Nx%, Ny%    : ' Coordonnees d'un point (pixels)
dim xt, yt      : ' Coordonnees d'un point (algebriques)
dim f_x, f_y    : ' Fonction complexe f(z)
dim fp_x, fp_y  : ' Derivee f'(z)

' --------------------------------------------------------------
' Variables globales de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

' Constantes mathematiques

dim MaxNum, MinNum, MaxLog, MinLog, Pi, PiDiv2

MaxLog =  709.78     : ' Argument max. pour EXP
MinLog = -708.39     : ' Argument min. pour EXP
MaxNum = exp(MaxLog) : ' Nb reel max. ~ 2^1024
MinNum = exp(MinLog) : ' Nb reel min. ~ 2^(-1022)
Pi     = 4 * atn(1)
PiDiv2 = Pi / 2

' Resultats des calculs
' Partie reelle, partie imaginaire, module, argument, signe

dim r_x, r_y, r_mod, r_arg, r_sgn

' Code d'erreur
'  0 = pas d'erreur
' -1 = argument hors bornes
' -2 = singularite
' -3 = overflow
' -4 = underflow

dim ErrCode%

' --------------------------------------------------------------
' Description des objets
' --------------------------------------------------------------

' Fenetre principale

left    0, 50
top     0, 50
width   0, PicWidth% + 70
height  0, PicHeight% + 120
caption 0, "Fractale de Newton ... Veuillez patienter."

' Image

picture 1
left    1, 30
top     1, 40
width   1, PicWidth%
height  1, PicHeight%

2d_target_is 1

' --------------------------------------------------------------
' Programme principal
' --------------------------------------------------------------

ColorFact = 0.01 * ColorFact
AbsColor  = abs(ColorFact)
ScaleFact = 4 / (PicHeight% * ZoomFact)

for Ny% = 0 to PicHeight% - 1
 yt = y0 - ScaleFact * (Ny% - HalfPicHeight%)
 for Nx% = 0 to PicWidth% - 1
   xt = x0 + ScaleFact * (Nx% - HalfPicWidth%)
   Newton(xt, yt)
   2d_pen_color R%, G%, B%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

file_save 1, "newton.bmp"

caption 0, "Fractale de Newton ... Terminé."

end

' --------------------------------------------------------------
' Sous-programmes
' --------------------------------------------------------------

sub Func(x, y)
' Calcul de la fonction complexe f(z) = (f_x, f_y)
' et de sa derivee f'(z) = (fp_x, fp_y)
' Ici f(z) = z^p - 1 et f'(z) = p z^(p - 1)

 dim_local u_x, u_y     : ' u = z^(p - 1)

 CIntPower(x, y, p% - 1)
 
 u_x = r_x
 u_y = r_y
 
 CMul(u_x, u_y, x, y)

 f_x = r_x - 1
 f_y = r_y

 fp_x = p% * u_x
 fp_y = p% * u_y
end_sub

sub Newton(xt, yt)
' Iteration de la fonction complexe au point (xt, yt)

 dim_local Iter%, mf2, mq2, x, y, zn_x, zn_y, q_x, q_y

 zn_x = xt
 zn_y = yt

 Iter% = 0

 repeat
   x = zn_x
   y = zn_y

   Func(x, y)

   if fp_x = 0 and fp_y = 0
     Iter% = MaxIter%             : ' Arret du calcul si f'(z) = 0
   else
     CDiv(f_x, f_y, fp_x, fp_y)   : ' q = f(z) / f'(z)
     
     q_x = r_x
     q_y = r_y

     zn_x = x - q_x
     zn_y = y - q_y

     mf2 = f_x * f_x + f_y * f_y  : ' |f(z)|^2
     mq2 = q_x * q_x + q_y * q_y  : ' |q|^2

     Iter% = Iter% + 1
   end_if
 until (mf2 < Eps2 and mq2 < Eps2) or (Iter% > MaxIter%)

 if Iter% > MaxIter%  : ' Pas de convergence
   R% = 0
   G% = 0
   B% = 0
 else
   MdbCol(Iter%)
 end_if
end_sub

sub MdbCol(Iter%)
' Determine la teinte (Hue, H) et la Saturation (S)
' d'apres le nombre d'iterations

 dim_local Q, Angle, Radius, H, S, V

 Q = log(Iter%) * AbsColor

 if Q < 0.5
   Q = 1 - 1.5 * Q
   Angle = 1 - Q
 else
   Q = 1.5 * Q - 0.5
   Angle = Q
 end_if

 Radius = sqr(Q)

' Si ColorFact > 0, assombrir une bande sur 2

 V = CstVal

 if (ColorFact > 0) and (odd(Iter%) > 0)
   V = 0.85 * V
   Radius = 0.667 * Radius
 end_if

 H = Angle * 10
 H = H - int(H)
 H = H * 360

 S = Radius - int(Radius)

' Convertir HSV en RGB

 HSVtoRGB(H, S, V)
end_sub

sub HSVtoRGB(H, S, V)
' Conversion HSV --> RGB. Resultats dans R%, G%, B%

 dim_local II%, ZZ, FF, PP, QQ, TT, RR, GG, BB

 if S = 0
   R% = int(V * 255)
   G% = R%
   B% = R%
 else
   ZZ  = H / 60
   II% = int(ZZ)
   FF  = ZZ - int(ZZ)
   PP  = V * (1 - S)
   QQ  = V * (1 - S * FF)
   TT  = V * (1 - S * (1 - FF))

   select II%
     case 0
       RR = V  : GG = TT : BB = PP
     case 1
       RR = QQ : GG = V  : BB = PP
     case 2
       RR = PP : GG = V  : BB = TT
     case 3
       RR = PP : GG = QQ : BB = V
     case 4
       RR = TT : GG = PP : BB = V
     case 5
       RR = V  : GG = PP : BB = QQ
   end_select

   R% = int(RR * 255)
   G% = int(GG * 255)
   B% = int(BB * 255)
 end_if
end_sub

' --------------------------------------------------------------
' Procedures de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

sub CMul(a_x, a_y, b_x, b_y)
' Multiplication : r_x + i r_y = (a_x + i a_y) * (b_x + i b_y)

 ErrCode% = 0

 r_x = a_x * b_x - a_y * b_y
 r_y = a_x * b_y + a_y * b_x
end_sub

sub CDiv(a_x, a_y, b_x, b_y)
' Division : r_x + i r_y = (a_x + i a_y) / (b_x + i b_y)
' Algorithme d'apres "Numerical Recipes"

 dim_local q, t

 if b_x = 0 and b_y = 0
   ErrCode% = -3
   r_x = MaxNum
   r_y = MaxNum
 else
   ErrCode% = 0
   if abs(b_x) >= abs(b_y)
     q = b_y / b_x
     t = b_x + b_y * q
     r_x = (a_x + a_y * q) / t
     r_y = (a_y - a_x * q) / t
   else
     q = b_x / b_y
     t = b_x * q + b_y
     r_x = (a_x * q + a_y) / t
     r_y = (a_y * q - a_x) / t
   end_if
 end_if
end_sub

sub CSquare(a_x, a_y)
' Carre : r_x + i r_y = (a_x + i a_y)^2

 ErrCode% = 0

 r_x = a_x * a_x - a_y * a_y
 r_y = 2 * a_x * a_y
end_sub

sub CIntPower(a_x, a_y, n%)
' Puissance entiere : r_x + i r_y = (a_x + i a_y)^n

 dim_local m%, b_x, b_y, res_x, res_y
 
 ErrCode% = 0
 
 if a_x = 0 and a_y = 0
   if n% = 0
     ' 0^0 = lim x^x quand x --> 0 = 1
     r_x = 1
     r_y = 0
   else
     if n% > 0
       ' 0^n = 0 si n > 0
       r_x = 0
       r_y = 0
     else
       ' 0^n indefini si n < 0
       ErrCode% = -2
       r_x = MaxNum
       r_y = MaxNum
     end_if
   end_if
 else
   if n% < 0
     m% = abs(n%)
     CInv(a_x, a_y)
     b_x = r_x
     b_y = r_y
   else
     m% = n%
     b_x = a_x
     b_y = a_y
   end_if

   res_x = 1 : res_y = 0

   while m% > 0
     if odd(m%) = 1
       CMul(b_x, b_y, res_x, res_y)
       res_x = r_x
       res_y = r_y
     end_if
     CSquare(b_x, b_y)
     b_x = r_x
     b_y = r_y
     m% = int(m% / 2)
   end_while
 
   r_x = res_x
   r_y = res_y
 end_if  
end_sub

' **************************************************************
' Fractale de Newton : Fonction z^p - 1 (p reel)
' **************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image
dim PicHeight% : PicHeight% = 480   : '   en pixels
dim p          : p = 3.5            : ' Exposant (z^p - 1)
dim x0         : x0 = 0             : ' Coordonnees du centre
dim y0         : y0 = 0             : '   de l'image
dim Eps        : Eps = 0.000001     : ' Precision requise
dim MaxIter%   : MaxIter% = 100     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 1       : ' Facteur de zoom
dim ColorFact  : ColorFact = -2     : ' Facteur de coloration
dim CstVal     : CstVal = 0.9       : ' Luminosite HSV constante

' Variables supplementaires

dim HalfPicWidth%  : HalfPicWidth%  = PicWidth% / 2
dim HalfPicHeight% : HalfPicHeight% = PicHeight% / 2
dim Eps2           : Eps2 = Eps * Eps

dim R%, G%, B%  : ' Parametres de coloration
dim AbsColor    : ' abs(ColorFact)
dim ScaleFact   : ' Facteur d'echelle = distance entre 2 pixels
dim Nx%, Ny%    : ' Coordonnees d'un point (pixels)
dim xt, yt      : ' Coordonnees d'un point (algebriques)
dim f_x, f_y    : ' Fonction complexe f(z)
dim fp_x, fp_y  : ' Derivee f'(z)

' --------------------------------------------------------------
' Variables globales de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

' Constantes mathematiques

dim MaxNum, MinNum, MaxLog, MinLog, Pi, PiDiv2

MaxLog =  709.78     : ' Argument max. pour EXP
MinLog = -708.39     : ' Argument min. pour EXP
MaxNum = exp(MaxLog) : ' Nb reel max. ~ 2^1024
MinNum = exp(MinLog) : ' Nb reel min. ~ 2^(-1022)
Pi     = 4 * atn(1)
PiDiv2 = Pi / 2

' Resultats des calculs
' Partie reelle, partie imaginaire, module, argument, signe

dim r_x, r_y, r_mod, r_arg, r_sgn

' Code d'erreur
'  0 = pas d'erreur
' -1 = argument hors bornes
' -2 = singularite
' -3 = overflow
' -4 = underflow

dim ErrCode%

' --------------------------------------------------------------
' Description des objets
' --------------------------------------------------------------

' Fenetre principale

left    0, 50
top     0, 50
width   0, PicWidth% + 70
height  0, PicHeight% + 120
caption 0, "Fractale de Newton ... Veuillez patienter."

' Image

picture 1
left    1, 30
top     1, 40
width   1, PicWidth%
height  1, PicHeight%

2d_target_is 1

' --------------------------------------------------------------
' Programme principal
' --------------------------------------------------------------

ColorFact = 0.01 * ColorFact
AbsColor  = abs(ColorFact)
ScaleFact = 4 / (PicHeight% * ZoomFact)

for Ny% = 0 to PicHeight% - 1
 yt = y0 - ScaleFact * (Ny% - HalfPicHeight%)
 for Nx% = 0 to PicWidth% - 1
   xt = x0 + ScaleFact * (Nx% - HalfPicWidth%)
   Newton(xt, yt)
   2d_pen_color R%, G%, B%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

file_save 1, "newton.bmp"

caption 0, "Fractale de Newton ... Terminé."

end

' --------------------------------------------------------------
' Sous-programmes
' --------------------------------------------------------------

sub Func(x, y)
' Calcul de la fonction complexe f(z) = (f_x, f_y)
' et de sa derivee f'(z) = (fp_x, fp_y)
' Ici f(z) = z^p - 1 et f'(z) = p z^(p - 1)

 dim_local u_x, u_y     : ' u = z^(p - 1)

 CRealPower(x, y, p - 1)
 
 u_x = r_x
 u_y = r_y
 
 CMul(u_x, u_y, x, y)

 f_x = r_x - 1
 f_y = r_y

 fp_x = p * u_x
 fp_y = p * u_y
end_sub

sub Newton(xt, yt)
' Iteration de la fonction complexe au point (xt, yt)

 dim_local Iter%, mf2, mq2, x, y, zn_x, zn_y, q_x, q_y

 zn_x = xt
 zn_y = yt

 Iter% = 0

 repeat
   x = zn_x
   y = zn_y

   Func(x, y)

   if fp_x = 0 and fp_y = 0
     Iter% = MaxIter%             : ' Arret du calcul si f'(z) = 0
   else
     CDiv(f_x, f_y, fp_x, fp_y)   : ' q = f(z) / f'(z)
     
     q_x = r_x
     q_y = r_y

     zn_x = x - q_x
     zn_y = y - q_y

     mf2 = f_x * f_x + f_y * f_y  : ' |f(z)|^2
     mq2 = q_x * q_x + q_y * q_y  : ' |q|^2

     Iter% = Iter% + 1
   end_if
 until (mf2 < Eps2 and mq2 < Eps2) or (Iter% > MaxIter%)

 if Iter% > MaxIter%  : ' Pas de convergence
   R% = 0
   G% = 0
   B% = 0
 else
   MdbCol(Iter%)
 end_if
end_sub

sub MdbCol(Iter%)
' Determine la teinte (Hue, H) et la Saturation (S)
' d'apres le nombre d'iterations

 dim_local Q, Angle, Radius, H, S, V

 Q = log(Iter%) * AbsColor

 if Q < 0.5
   Q = 1 - 1.5 * Q
   Angle = 1 - Q
 else
   Q = 1.5 * Q - 0.5
   Angle = Q
 end_if

 Radius = sqr(Q)

' Si ColorFact > 0, assombrir une bande sur 2

 V = CstVal

 if (ColorFact > 0) and (odd(Iter%) > 0)
   V = 0.85 * V
   Radius = 0.667 * Radius
 end_if

 H = Angle * 10
 H = H - int(H)
 H = H * 360

 S = Radius - int(Radius)

' Convertir HSV en RGB

 HSVtoRGB(H, S, V)
end_sub

sub HSVtoRGB(H, S, V)
' Conversion HSV --> RGB. Resultats dans R%, G%, B%

 dim_local II%, ZZ, FF, PP, QQ, TT, RR, GG, BB

 if S = 0
   R% = int(V * 255)
   G% = R%
   B% = R%
 else
   ZZ  = H / 60
   II% = int(ZZ)
   FF  = ZZ - int(ZZ)
   PP  = V * (1 - S)
   QQ  = V * (1 - S * FF)
   TT  = V * (1 - S * (1 - FF))

   select II%
     case 0
       RR = V  : GG = TT : BB = PP
     case 1
       RR = QQ : GG = V  : BB = PP
     case 2
       RR = PP : GG = V  : BB = TT
     case 3
       RR = PP : GG = QQ : BB = V
     case 4
       RR = TT : GG = PP : BB = V
     case 5
       RR = V  : GG = PP : BB = QQ
   end_select

   R% = int(RR * 255)
   G% = int(GG * 255)
   B% = int(BB * 255)
 end_if
end_sub

' --------------------------------------------------------------
' Procedures de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

sub CMul(a_x, a_y, b_x, b_y)
' Multiplication : r_x + i r_y = (a_x + i a_y) * (b_x + i b_y)

 ErrCode% = 0

 r_x = a_x * b_x - a_y * b_y
 r_y = a_x * b_y + a_y * b_x
end_sub

sub CDiv(a_x, a_y, b_x, b_y)
' Division : r_x + i r_y = (a_x + i a_y) / (b_x + i b_y)
' Algorithme d'apres "Numerical Recipes"

 dim_local q, t

 if b_x = 0 and b_y = 0
   ErrCode% = -3
   r_x = MaxNum
   r_y = MaxNum
 else
   ErrCode% = 0
   if abs(b_x) >= abs(b_y)
     q = b_y / b_x
     t = b_x + b_y * q
     r_x = (a_x + a_y * q) / t
     r_y = (a_y - a_x * q) / t
   else
     q = b_x / b_y
     t = b_x * q + b_y
     r_x = (a_x * q + a_y) / t
     r_y = (a_y * q - a_x) / t  
   end_if    
 end_if
end_sub

sub CAbs(a_x, a_y)
' Module : r_mod = |a_x + i a_y|
' Algorithme d'apres "Numerical Recipes"

 ErrCode% = 0

 dim_local AbsX, AbsY, R, C

 AbsX = abs(a_x)
 AbsY = abs(a_y)

 if a_x = 0
   r_mod = abs(a_y)
 else
   if a_y = 0
     r_mod = abs(a_x)
   else
     if AbsX > AbsY
       R = AbsY / AbsX
       C = AbsX
     else
       R = AbsX / AbsY
       C = AbsY
     end_if
     r_mod = C * sqr(1 + R * R)
   end_if
 end_if
end_sub

sub CArg(a_x, a_y)
' Argument : r_arg = arg(a_x + i a_y)
' Resultat dans [-Pi, Pi)
' Equivaut a atan2(a_y, a_x)

 ErrCode% = 0

 if a_x = 0
   r_arg = sgn(a_y) * PiDiv2
 else
   ' 4e / 1er quadrant : -Pi/2..Pi/2
   r_arg = atn(a_y / a_x)
   if a_x < 0
     if a_y > 0
       ' 2e quadrant : Pi/2..Pi
       r_arg = r_arg + Pi
     else
       ' 3e quadrant : -Pi..-Pi/2
       r_arg = r_arg - Pi
     end_if
   end_if
 end_if
end_sub

sub CRealPower(a_x, a_y, p)
' Puissance (exposant reel) : (a_x + i a_y)^p
' Resultat dans r_x, r_y
' Resultat aussi dans r_mod, r_arg si a <> 0

 ErrCode% = 0

 if a_x = 0 and a_y = 0
   if p = 0
     ' 0^0 = lim x^x quand x --> 0 = 1
     r_x = 1
     r_y = 0
   else
     if p > 0
       ' 0^p = 0 si p > 0
       r_x = 0
       r_y = 0
     else
       ' 0^p indefini si p < 0
       ErrCode% = -2
       r_x = MaxNum
       r_y = MaxNum
     end_if
   end_if
 else
   CAbs(a_x, a_y)
   CArg(a_x, a_y)
   r_mod = power(r_mod, p)
   r_arg = r_arg * p
   r_x = r_mod * cos(r_arg)
   r_y = r_mod * sin(r_arg)
 end_if
end_sub

' **************************************************************
' Fractale de Newton : Fonction z^p - 1 (p complexe)
' **************************************************************

' Variables definies par l'utilisateur

dim PicWidth%  : PicWidth%  = 640   : ' Taille de l'image
dim PicHeight% : PicHeight% = 480   : '   en pixels
dim p_x, p_y   : p_x = 3 : p_y = 2  : ' Exposant (z^p - 1)
dim x0         : x0 = 0             : ' Coordonnees du centre
dim y0         : y0 = 0             : '   de l'image
dim Eps        : Eps = 0.000001     : ' Precision requise
dim MaxIter%   : MaxIter% = 100     : ' Nb maxi d'iterations
dim ZoomFact   : ZoomFact = 1       : ' Facteur de zoom
dim ColorFact  : ColorFact = -2     : ' Facteur de coloration
dim CstVal     : CstVal = 0.9       : ' Luminosite HSV constante

' Variables supplementaires

dim HalfPicWidth%  : HalfPicWidth%  = PicWidth% / 2
dim HalfPicHeight% : HalfPicHeight% = PicHeight% / 2
dim Eps2           : Eps2 = Eps * Eps

dim R%, G%, B%  : ' Parametres de coloration
dim AbsColor    : ' abs(ColorFact)
dim ScaleFact   : ' Facteur d'echelle = distance entre 2 pixels
dim Nx%, Ny%    : ' Coordonnees d'un point (pixels)
dim xt, yt      : ' Coordonnees d'un point (algebriques)
dim f_x, f_y    : ' Fonction complexe f(z)
dim fp_x, fp_y  : ' Derivee f'(z)

' --------------------------------------------------------------
' Variables globales de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

' Constantes mathematiques

dim MaxNum, MinNum, MaxLog, MinLog, Pi, PiDiv2

MaxLog =  709.78     : ' Argument max. pour EXP
MinLog = -708.39     : ' Argument min. pour EXP
MaxNum = exp(MaxLog) : ' Nb reel max. ~ 2^1024
MinNum = exp(MinLog) : ' Nb reel min. ~ 2^(-1022)
Pi     = 4 * atn(1)
PiDiv2 = Pi / 2

' Resultats des calculs
' Partie reelle, partie imaginaire, module, argument, signe

dim r_x, r_y, r_mod, r_arg, r_sgn

' Code d'erreur
'  0 = pas d'erreur
' -1 = argument hors bornes
' -2 = singularite
' -3 = overflow
' -4 = underflow

dim ErrCode%

' --------------------------------------------------------------
' Description des objets
' --------------------------------------------------------------

' Fenetre principale

left    0, 50
top     0, 50
width   0, PicWidth% + 70
height  0, PicHeight% + 120
caption 0, "Fractale de Newton ... Veuillez patienter."

' Image

picture 1
left    1, 30
top     1, 40
width   1, PicWidth%
height  1, PicHeight%

2d_target_is 1

' --------------------------------------------------------------
' Programme principal
' --------------------------------------------------------------

ColorFact = 0.01 * ColorFact
AbsColor  = abs(ColorFact)
ScaleFact = 4 / (PicHeight% * ZoomFact)

for Ny% = 0 to PicHeight% - 1
 yt = y0 - ScaleFact * (Ny% - HalfPicHeight%)
 for Nx% = 0 to PicWidth% - 1
   xt = x0 + ScaleFact * (Nx% - HalfPicWidth%)
   Newton(xt, yt)
   2d_pen_color R%, G%, B%
   2d_point Nx%, Ny%
 next Nx%
next Ny%

file_save 1, "newton.bmp"

caption 0, "Fractale de Newton ... Terminé."

end

' --------------------------------------------------------------
' Sous-programmes
' --------------------------------------------------------------

sub Func(x, y)
' Calcul de la fonction complexe f(z) = (f_x, f_y)
' et de sa derivee f'(z) = (fp_x, fp_y)
' Ici f(z) = z^p - 1 et f'(z) = p z^(p - 1)

 dim_local u_x, u_y        : ' u = z^(p - 1)

 CPower(x, y, p_x - 1, p_y)
 
 u_x = r_x
 u_y = r_y
 
 CMul(u_x, u_y, x, y)      : ' z^p

 f_x = r_x - 1
 f_y = r_y
 
 CMul(p_x, p_y, u_x, u_y)  : ' p z^(p - 1)

 fp_x = r_x
 fp_y = r_y
end_sub

sub Newton(xt, yt)
' Iteration de la fonction complexe au point (xt, yt)

 dim_local Iter%, mf2, mq2, x, y, zn_x, zn_y, q_x, q_y

 zn_x = xt
 zn_y = yt

 Iter% = 0

 repeat
   x = zn_x
   y = zn_y

   Func(x, y)

   if fp_x = 0 and fp_y = 0
     Iter% = MaxIter%             : ' Arret du calcul si f'(z) = 0
   else
     CDiv(f_x, f_y, fp_x, fp_y)   : ' q = f(z) / f'(z)
     
     q_x = r_x
     q_y = r_y

     zn_x = x - q_x
     zn_y = y - q_y

     mf2 = f_x * f_x + f_y * f_y  : ' |f(z)|^2
     mq2 = q_x * q_x + q_y * q_y  : ' |q|^2

     Iter% = Iter% + 1
   end_if
 until (mf2 < Eps2 and mq2 < Eps2) or (Iter% > MaxIter%)

 if Iter% > MaxIter%  : ' Pas de convergence
   R% = 0
   G% = 0
   B% = 0
 else
   MdbCol(Iter%)
 end_if
end_sub

sub MdbCol(Iter%)
' Determine la teinte (Hue, H) et la Saturation (S)
' d'apres le nombre d'iterations

 dim_local Q, Angle, Radius, H, S, V

 Q = log(Iter%) * AbsColor

 if Q < 0.5
   Q = 1 - 1.5 * Q
   Angle = 1 - Q
 else
   Q = 1.5 * Q - 0.5
   Angle = Q
 end_if

 Radius = sqr(Q)

' Si ColorFact > 0, assombrir une bande sur 2

 V = CstVal

 if (ColorFact > 0) and (odd(Iter%) > 0)
   V = 0.85 * V
   Radius = 0.667 * Radius
 end_if

 H = Angle * 10
 H = H - int(H)
 H = H * 360

 S = Radius - int(Radius)

' Convertir HSV en RGB

 HSVtoRGB(H, S, V)
end_sub

sub HSVtoRGB(H, S, V)
' Conversion HSV --> RGB. Resultats dans R%, G%, B%

 dim_local II%, ZZ, FF, PP, QQ, TT, RR, GG, BB

 if S = 0
   R% = int(V * 255)
   G% = R%
   B% = R%
 else
   ZZ  = H / 60
   II% = int(ZZ)
   FF  = ZZ - int(ZZ)
   PP  = V * (1 - S)
   QQ  = V * (1 - S * FF)
   TT  = V * (1 - S * (1 - FF))

   select II%
     case 0
       RR = V  : GG = TT : BB = PP
     case 1
       RR = QQ : GG = V  : BB = PP
     case 2
       RR = PP : GG = V  : BB = TT
     case 3
       RR = PP : GG = QQ : BB = V
     case 4
       RR = TT : GG = PP : BB = V
     case 5
       RR = V  : GG = PP : BB = QQ
   end_select

   R% = int(RR * 255)
   G% = int(GG * 255)
   B% = int(BB * 255)
 end_if
end_sub

' --------------------------------------------------------------
' Procedures de la bibliotheque (nombres complexes)
' --------------------------------------------------------------

sub CMul(a_x, a_y, b_x, b_y)
' Multiplication : r_x + i r_y = (a_x + i a_y) * (b_x + i b_y)

 ErrCode% = 0

 r_x = a_x * b_x - a_y * b_y
 r_y = a_x * b_y + a_y * b_x
end_sub

sub CDiv(a_x, a_y, b_x, b_y)
' Division : r_x + i r_y = (a_x + i a_y) / (b_x + i b_y)
' Algorithme d'apres "Numerical Recipes"

 dim_local q, t

 if b_x = 0 and b_y = 0
   ErrCode% = -3
   r_x = MaxNum
   r_y = MaxNum
 else
   ErrCode% = 0
   if abs(b_x) >= abs(b_y)
     q = b_y / b_x
     t = b_x + b_y * q
     r_x = (a_x + a_y * q) / t
     r_y = (a_y - a_x * q) / t
   else
     q = b_x / b_y
     t = b_x * q + b_y
     r_x = (a_x * q + a_y) / t
     r_y = (a_y * q - a_x) / t  
   end_if    
 end_if
end_sub

sub CAbs(a_x, a_y)
' Module : r_mod = |a_x + i a_y|
' Algorithme d'apres "Numerical Recipes"

 ErrCode% = 0

 dim_local AbsX, AbsY, R, C

 AbsX = abs(a_x)
 AbsY = abs(a_y)

 if a_x = 0
   r_mod = abs(a_y)
 else
   if a_y = 0
     r_mod = abs(a_x)
   else
     if AbsX > AbsY
       R = AbsY / AbsX
       C = AbsX
     else
       R = AbsX / AbsY
       C = AbsY
     end_if
     r_mod = C * sqr(1 + R * R)
   end_if
 end_if
end_sub

sub CArg(a_x, a_y)
' Argument : r_arg = arg(a_x + i a_y)
' Resultat dans [-Pi, Pi)
' Equivaut a atan2(a_y, a_x)

 ErrCode% = 0

 if a_x = 0
   r_arg = sgn(a_y) * PiDiv2
 else
   ' 4e / 1er quadrant : -Pi/2..Pi/2
   r_arg = atn(a_y / a_x)
   if a_x < 0
     if a_y > 0
       ' 2e quadrant : Pi/2..Pi
       r_arg = r_arg + Pi
     else
       ' 3e quadrant : -Pi..-Pi/2
       r_arg = r_arg - Pi
     end_if
   end_if
 end_if
end_sub

sub CExp(a_x, a_y)
' Exponentielle complexe : r_x + i r_y = exp(a_x + i a_y)

 dim_local ExpX

 if a_x < MinLog
   ErrCode% = -4
   r_x = 0
   r_y = 0
 else
   if a_x > MaxLog
     ErrCode = -3
     ExpX = MaxNum
   else
     ErrCode% = 0
     ExpX = exp(a_x)
   end_if
   r_x = ExpX * cos(a_y)
   r_y = ExpX * sin(a_y)
 end_if
end_sub

sub CPower(a_x, a_y, b_x, b_y)
' Puissance (exposant complexe) : (a_x + i a_y)^(b_x + i b_y)
' Resultat dans r_x, r_y

 ErrCode% = 0

 if a_x = 0 and a_y = 0
   if b_x = 0 and b_y = 0
     ' 0^0 = lim x^x quand x --> 0 = 1
     r_x = 1
     r_y = 0
   else
     ' 0^p = 0 si p > 0
     r_x = 0
     r_y = 0
   end_if
 else
   ' exp(b ln(a))
   CAbs(a_x, a_y)
   CArg(a_x, a_y)
   CMul(b_x, b_y, log(r_mod), r_arg)
   CExp(r_x, r_y)
 end_if
end_sub