As I use the terms, the dual numbers are a different thing from non-standard analysis. Non-standard analysis, as I understand the term, uses non-standard models of the real numbers, and its infinitesimals do not satisfy h^2 = 0. In non-standard analysis, f'(x) is the standard part of (f(x+h)-f(x))/h , for an infinitesimal h (i.e. for a non-standard real which is smaller than any non-zero standard rational number). (In order to apply this definition, f should be defined in a way which does not use anything requiring determining if a number is standard, or taking the standard part of something, etc.)

The dual numbers, on the other hand, are the ring R[h]/(h^2) . This is not a field, while non-standard models of the real numbers do form fields.

The dual numbers suffice to define differentiation of polynomials (which may not be sufficient for some purposes! [a]), and something like dual numbers is used in algebraic geometry to define the Zariski tangent spaces for points of algebraic varieties (whether in characteristic 0 or in characteristic p. In characteristic p, one certainly can't use an epsilon-delta definition!).

I really don't see your point about the "gee-whiz factor". While different things, both non-standard analysis and dual numbers can be handled rigorously, and have their use-cases, even though I certainly would at least default to thinking of differentiation in terms of the limits definition (assuming I'm thinking of any specific definition at all).

I assume that the counter-example you refer to is the (dx/dy)(dy/dz)(dz/dx) thing. That indeed doesn't seem like the kind of thing that using the dual numbers would be especially suited for. Though, also, not the sort of thing that should really come up in auto-diff I would think?

If there is a common parameterization of the values of x,y,z by some variable t (on some interval I), on some neighborhood of the point under consideration, where {(x(t),y(t)) | t in I}, {(y(t),z(t)) | t in I}, and {(z(t),x(t)) | t in I}, are each differentiable functions, and where they satisfy the relationships between the variables x,y,z from the larger context, and all three of (dx/dy), (dy/dz), (dz/dx) exist at the point in question, then it seems that the product should be 1.

[a] Though, if the function is analytic (not just smooth), unless I'm missing something, it should also give the right answer (but still not a good definition of course, because should define differentiation before defining what it is for a function to be analytic.) (of course, just because the function is analytic doesn't make its domain include the dual numbers. One has to take a power series for it, and apply this power series to the element of the ring of dual numbers, not apply the original function to it.)