Help for generating Jenkins-Traub in Ultra Fractal

[avachon100510]

I could really use some of your expertise to help me make this work here:

Jenkins-Traub_CubeRootsofUnity {
init:
  import "common.ulb"
  z = 0.5
  s = #pixel
loop:


  complex func p(complex x)
    return x^3 - 1
  endfunc
  complex func dp(complex x)
    return 3*x^2
  endfunc

  complex func update_H(complex H_lambda_x, complex H_lambda_s, complex x, complex s)
    complex p_x = p(x)
    complex p_s = p(s)
    complex H_lambda_x_new = p_x - (p_s/H_lambda_s)*H_lambda_x
    complex H_lambda_x = 1/(x - s)*H_lambda_x_new
    return H_lambda_x
  endfunc

  int m = @num_maxlambda_Stage1
  ComplexArray H = new ComplexArray(m)
  ComplexArray H0 = new ComplexArray(m)

  H.m_Elements[0] = dp(z)
  H0.m_Elements[0] = dp(s)
  int lambda = 0

 while (lambda <= m)
   H.m_Elements[lambda + 1] = update_H(H.m_Elements[lambda], H0.m_Elements[lambda], z, s)
   olds = s
   s = s - p(s)/(H0.m_Elements[lambda])
   H0.m_Elements[lambda + 1] = H.m_Elements[lambda + 1]
   lambda = lambda + 1
 endwhile
  oldz = z
  z = z - p(z)/(H.m_Elements[lambda])
bailout:
   |s - olds| >= @bailout && |p(s)| >= @bailout
default:
  title = "Jenkins-Traub, cube roots of unity"
param bailout
  caption = "Bailout value"
  default = 0.000001
endparam
int param num_maxlambda_Stage1
  caption = "Maximum G_lambda iterations"
  default = 5
  hint = "Usually, m = 5 is chosen for polynomials of degree up to 50. Only change this if you use polynomials of higher degree."
endparam

}

[osveliz]

I haven't played with UltraFractal very much but I don't know how well it will do with multiplying and diving polynomials. For example in my function call shift in the example code:

def shift(P: Poly, H: Poly, s=0):
    """
    Shifts a polynomial H by evaluating it at a shift point s using polynomial division.

    Parameters
    ----------
    P : Poly
        The input polynomial.
    H : Poly
        The polynomial to be shifted.
    s : float, optional
        The shift point (default is 0).

    Returns
    -------
    Poly
        The shifted H.
    """
    Ps = P.eval(s).evalf()
    Hs = H.eval(s).evalf()
    result, _ = div(H-P*Hs/Ps, x-s)
    return result

The variables P and H are SymPy Polynomials so you can do things like evaluating the polynomial at a point Ps = P.eval(s) and multiplying polynomials agains scalars like P*Hs but I also needed to do polynomial division like in the second to last line result, _ = div(H-P*Hs/Ps, x-s) simply applying Horner's Method to divide the resulting H-P*Hs/Ps polynomial by x-2 and dropping the remainder.

I also don't know if writing something in UltraFractal like complex p_x = p(x) does what I think you intend for it to do. I suspect that the complex type only allows you to store a number, not a function.

I might suggest writing your code to just use arrays of coefficients that you can evaluate with Horner and also to do things like multiplying by a scalar (multiply every coefficient by that value) and also allow for subtracting polynomials (subtract by index and account for if the indexes are different sizes).

[avachon100510]

Arrays are something I did consider at one point. But the thing is: I want to generate a fractal as well and try to experiment with it. I'd need to craft a Polynomial class or something like a function class really. The parent function would have to store an input variable, and be updated without affecting the input variable(s) themselves, and part of the H polynomial updating process.

[osveliz]

Wouldn't a Polynomial class (or library of functions / something like it) still essentially implement a polynomial via an array anyway? I don't think there's another good way to save an unknown number of coefficients for any generic algebraic equation.

[avachon100510]

Well, this is why the Array classes in the common.ulb plug-in library exist. Since arrays themselves cannot be taken as a function argument, they come in handy!