Hi @vegardjervell - could you please try the following:

dual z_func(dual a){
    const auto rootfun = [&](dual z_) -> dual {return pow(z_, 2.) - a;};
    const auto rootfun_z = [&](dual z_) -> dual {return 2*z_;};
    dual f0, f1;
    dual z = 1;
    do {
        f0 = rootfun(z);
        f1 = rootfun_z(z);  // instead of derivative(rootfun, wrt(z), at(z)), which breaks derivative propagation because the result is not dual but an order lower (double)
        z -= f0 / f1;
    } while (abs(f0) > 1e-6);
    return z;
}

int main(){
    dual a = 4.;
    auto [z0, z1] = derivatives(z_func, wrt(a), at(a));
    std::cout << "Z : " << z0 << ", " << z1 << std::endl; // Should be (2, 0.25) 
    return 0;
}

Hopefully this compiles - cannot check now.

Explanation in the comment above. You are "consuming" the autodiff::dual number in the loop, producing a double, and this breaks the chain of derivative propagation.

Thank you, your example compiles and works, but has the issue that it requires differentiating rootfun "manually". In this sample case that's trivial, but in my actual use-case I would very much like to avoid having to do that. Taking your advice regarding how the dual is consumed, I ended up with this working example:

dual z_func(dual2nd a){
    const auto rootfun = [&](dual2nd z_) -> dual2nd {return pow(z_, 2.) - a;};
    const auto rootfun_z = [&](dual2nd z_) -> dual {return derivative(rootfun, wrt(z_), at(z_));};
    dual2nd f0;
    dual f1;
    dual2nd z = 1;
    do {
        f0 = rootfun(z);
        f1 = rootfun_z(z);
        z -= f0 / f1;
    } while (abs(f0) > 1e-6);
    return z.val;
}

but ended up with some follow-up questions:

1. If I change f1 = rootfun_z(z) to f1 = derivative(rootfun, wrt(z), at(z)) (i.e. drop the lambda), this doesn't work. The same is true if I drop the explicit -> dual in the declaration of the rootfun_z lambda. I suspect there's some implicit conversion going on here that prevents the derivative from propagating properly, but can't tell immediately what it is.
2. If I change the initial guess for the solver to z = - 1, the solver converges correctly to $z = -2$, but an incorrect derivative (Analytically, $z^2 - a = 0 \implies \frac{dz}{da} = \frac{1}{2z}$, but at the solution $z = - 2$ I get the derivative -0.2667) I've tried reducing the tolerance and playing around with the initial guess, and found that the derivative will become incorrect if there are too few iterations (i.e. setting an initial guess far away from the solution leads to a correct derivative, while using an initial guess close to the solution leads to an incorrect derivative).
3. Perhaps most confusingly: If I so much as call derivative(rootfun, wrt(z), at(z)) in the loop (even without using it) I get the wrong derivative:

dual z_func(dual2 a){
    const auto rootfun = [&](dual2 z_) -> dual2 {return pow(z_, 2.) - a;};
    const auto rootfun_z = [&](dual2 z_) -> dual {return derivative(rootfun, wrt(z_), at(z_));};
    dual2 f0;
    dual f1;
    dual2 z = 1;
    do {
        f0 = rootfun(z);
        f1 = rootfun_z(z);
        dual unused = derivative(rootfun, wrt(z), at(z)); // Breaks expected behaviour, remove this line to get correct derivative of z_func
        std::cout << unused << std::endl; // ensure that "unused" is not optimized out.
        z -= f0 / f1;
    } while (abs(f0) > 1e-6);
    return z.val;
}

It definitely feels like there's a lot going on under the hood here that I need to account for, but I'm having a really hard time figuring out what the system to it is. Any help would be greatly appreciated!

I've done some more digging @allanleal, and I'm starting to wonder if this should be raised as an issue.

I'm not sure about the details of how this is implemented, but by comparing to my own (very bare-bones) implementation of hyperdual numbers it appears that the different "parts" of the dual numbers are being mixed up when I differentiate a function that needs to compute derivatives internally like this.

To give an example:

dual z_func(dual2nd a){
    const auto rootfun = [&](const dual2& z_, const dual2& a_) -> dual2 {return pow(z_, 2.) - a_;};
    dual2nd f0;
    dual f1;
    dual2nd z = 1;
    dual2nd true_dfdz;
    do {
        f0 = rootfun(z, a);
        f1 = derivative(rootfun, wrt(z), at(z, a));
        true_dfdz = 2. * z;
        std::cout << "z = " << z << ", f(z) = " << f0 << ", f'(z) = " << f1 << " / " << true_dfdz << std::endl;
        z -= f0 / f1;
    } while (abs(f0) > 1e-6);
    return z.val;
}

int main(){
    dual2nd a = 4.;
    dual z0 = z_func(a);
    std::cout << "-------------" << std::endl;
    double z1 = derivative(z_func, wrt(a), at(a));
    return 0;
}

Gives the output

z = 1, f(z) = -3, f'(z) = 2 / 2
z = 2.5, f(z) = 2.25, f'(z) = 5 / 5
z = 2.05, f(z) = 0.2025, f'(z) = 4.1 / 4.1
z = 2.00061, f(z) = 0.0024394, f'(z) = 4.00122 / 4.00122
z = 2, f(z) = 3.71689e-07, f'(z) = 4 / 4
-------------
z = 1, f(z) = -3, f'(z) = 1 / 2
z = 4, f(z) = 12, f'(z) = 7 / 8
z = 2.28571, f(z) = 1.22449, f'(z) = 3.57143 / 4.57143
z = 1.94286, f(z) = -0.225306, f'(z) = 2.88571 / 3.88571
z = 2.02093, f(z) = 0.0841723, f'(z) = 3.04187 / 4.04187
z = 1.99326, f(z) = -0.0269056, f'(z) = 2.98652 / 3.98652
z = 2.00227, f(z) = 0.00909015, f'(z) = 3.00454 / 4.00454
z = 1.99925, f(z) = -0.00301632, f'(z) = 2.99849 / 3.99849
z = 2.00025, f(z) = 0.00100696, f'(z) = 3.0005 / 4.0005
z = 1.99992, f(z) = -0.000335483, f'(z) = 2.99983 / 3.99983
z = 2.00003, f(z) = 0.000111846, f'(z) = 3.00006 / 4.00006
z = 1.99999, f(z) = -3.72801e-05, f'(z) = 2.99998 / 3.99998
z = 2, f(z) = 1.24269e-05, f'(z) = 3.00001 / 4.00001
z = 2, f(z) = -4.14228e-06, f'(z) = 3 / 4
z = 2, f(z) = 1.38076e-06, f'(z) = 3 / 4
z = 2, f(z) = -4.60254e-07, f'(z) = 3 / 4

As you can see, the derivative computed inside z_func is correct when calling z_func directly, but not when calling derivative(z_func, wrt(a), at(a)), in which case it is offset by 1 compared to the correct derivative. I've tried various ways of defining rootfun (writing it as a stand-alone function, only passing z as an explicit variable, passing z and / or a as value vs. reference) but everything seems to give the same result: The derivative computed internally in z_func comes out wrong if z_func is called using derivative.

What is your take on this? For the time being it looks to me like I'll have to resort to my own implementation of hyperduals for these use-cases, is this something that should be raised as an issue?