(28/48/50) Dual Number Functions

[Albert Chan]

(03-20-2024 07:14 PM)John Keith Wrote:  Ángel Martin’s Double ROM manual has programs for dual number AGM and E(k). The example on p.34 gives a value of -2.622842115+ε*2.858860658 for the input 0.5+ε, whereas with the same input I get the value 1.68575035482+ε*.541731848587 with both DNEK and DELK.

Hi, John Keith

Your numbers are correct.

lua> D = require'dual'.D
lua> k = D.new(0.5, 1)
lua> a, b = 1+k, 1-k
lua> repeat b0=b[1]; a,b = (a+b):mul(.5), (a*b):sqrt(); until b0==b[1]
lua> unpack(b) -- = AGM(a,b)
0.9318083916224482      -0.29944545531661104
lua> K = pi/2 / b
lua> unpack(K) -- = K( k=0.5+ε )
1.685750354812596        0.5417318486132805

With dual numbers, we don't need modified AGM. We can get E from K
From https://en.wikipedia.org/wiki/Elliptic_integral

(dE/dk) = (E-K)/k
(dK/dk) = E/(k*(1-k*k)) - K/k = (dE/dk) / (1-k*k) + K * k/(1-k*k)

(dE/dk) = (1-k*k)*(dK/dk) - k*K
E = K + k*(dE/dk)

lua> E = D.new(0, (1-k[1]*k[1])*K[2] - k[1]*K[1])
lua> E[1] = K[1] + k[1]*E[2]
lua> unpack(E) -- = E( k=0.5+ε )
1.4674622093394272      -0.43657629094633765