programming - How to implement dual numbers in Mathematica?

I wonder how can I implement dual numbers in Mathematica, so that all functions work well with them (as with complex numbers).

Particularly, for each function $f$, $f(\varepsilon)=f^\prime(0)\varepsilon$

Answer

If you use the matrix representation, addition and multiplication works by default, but you need to use matrix multiplication and not just space (which is element-wise multiplication), for instance f[a_,b_]:=a.b+b.b would for instance work exactly like you would expect if a and b where dual numbers.

If you just want to have the rule enforce that you mention, you could just write that rule out:

 dualE /: f_[dualE] := f'[0] dualE

You can add other rules to make your other algebra work out aswell for intance for multiplication:

 dualE /: Power[dualE, 2] := 0
 (a + dualE b) (c + dualE d) // Expand // Collect[#, dualE] &

Update

If you want the full function application rule, then there is a problem with just defining it, lets assume you wanted to write down: f_[a+b dualE]:=f[a]+b f'[0] dualE if we expand this into the expression form it reads: f_[Plus[a,Times[b,dualE]]]:=... now such a definition can't have it's pattern attached to dualE, since it doesn't appear in the top level. What you would then need is to create a symbol that represents a dual number with both parts included, as Sasha suggested in a comment. Here is an example writing up rules for Times and Plus:

 dualE:=dualNumber[0,1];

 Times[dualNumber[r1_,d1_],dualNumber[r2_,d2_]] ^:= dualNumber[r1 r2,d2 r1+d1 r2]
 Times[dualNumber[r_,d_],n_] ^:= dualNumber[n r,n d]
 Plus[dualNumber[r_,d_],n_] ^:= dualNumber[r+n,d]

Then you can write the definition for function application with respect to dualNumber

 dualNumber /: f_[dualNumber[r_, d_]] := f[r] + d f'[r] dualE

And if you would like the numbers to be printed like 2+ 3 duelE rather then as dualNumber[2,3]. Then you can define a MakeBoxes rule for it.

 MakeBoxes[dualNumber[a_, d_], StandardForm] ^:= 
     RowBox[{MakeBoxes@a, "+", MakeBoxes@(Times[d, "dualE"])}]

Now then this works together to allow:

 Sin[(4 + dualE) (2 + 3 dualE)]

Sin[8] + dualE 14

 Sin[(a + b dualE) (c + d dualE)]

Sin[a c] + (b c + a d) dualE

Edited to correct error caught by Rahul Narain, also to improve incorrect output formating in symbolic expressions.