(28/48/50) Dual Number Functions

[John Keith]

The programs listed here provide functions for calculations with dual numbers, which are similar to complex numbers but instead of the imaginary unit i they use the infinitesimal unit ε. The programs are mostly simple translations of functions in Thomas Klemm's HP42S program here and in Ángel Martin’s Double ROM.

These programs use complex numbers as a convenient stand-in for dual numbers but these numbers are not actually complex numbers, and mixing the two will most likely yield nonsense results. However, the operations of addition, subtraction and negation are the same as for complex numbers, and the absolute value of a dual number is just the real part. Thus there are no programs for these functions.

The first listing is a directory of programs providing arithmetic and transcendental functions, as well as the program DC→R which is a modification of C→R which accepts either real or complex numbers. This is to allow easy mixing of dual numbers and real numbers, which are interpreted as (x, 0).

Code:


@ Dual number functions.
DIR
  DMUL @ Multiply two dual numbers.
  \<< DC\->R ROT DC\->R
      OVER 5 PICK * 5 ROLLD
      ROT ROT * ROT ROT * + R\->C
  \>>
  DDIV @ Divide two dual numbers.
  \<< DC\->R ROT DC\->R
      4 PICK * OVER 4 ROLL * -
      SWAP 3 PICK / ROT ROT SWAP SQ / R\->C
  \>>
  DINV @ Inverse (reciprocal) of a dual number.
  \<< DC\->R OVER INV ROT ROT
      NEG SWAP SQ / R\->C
  \>>
  DSQ @ Square of a dual number.
  \<< DC\->R OVER SQ ROT ROT
      * 2 * R\->C
  \>>
  DSQRT @ Square root of a dual number.
  \<< DC\->R SWAP \v/ SWAP OVER 2 * / R\->C
  \>>
  DLN @ Natural logarithm of a dual number.
  \<< DC\->R OVER / SWAP LN SWAP R\->C
  \>>
  DEXP @ Exponentiation (natural antilog) of a dual number.
  \<< DC\->R SWAP EXP SWAP OVER * R\->C
  \>>
  DIPWR @ Dual number to integer power.
  \<< SWAP DC\->R OVER 4 PICK ^
      4 ROLLD SWAP 3 PICK 1 - ^ * * R\->C
  \>>
  DDPWR @ Dual number to dual number power.
  \<< DC\->R ROT DC\->R OVER 5 PICK ^ 5 ROLLD
      OVER 5 ROLL SWAP OVER 1 - ^ * *
      SWAP LN ROT * 3 PICK * + R\->C
  \>>
  DNROOT @ Nth root of dual number for integer n.
  \<< SWAP DC\->R OVER 4 PICK INV ^
      4 ROLLD SWAP ABS 3 PICK INV 1 - ^
      * SWAP / R\->C
  \>>
  DC\->R @ Safe C-\>R for either real or dual number.
  \<< DUP TYPE 1 SAME
    \<< C\->R
    \>> 0 IFTE
  \>>
END

The next listing is the same directory optimized for the HP 49/50 using newer stack commands and decimal points for real numbers.

Code:


@ Dual number functions.
DIR
  DMUL @ Multiply two dual numbers.
  \<< DC\->R ROT DC\->R
      OVER 5. PICK * 5. ROLLD
      UNROT * UNROT * + R\->C
  \>>
  DDIV @ Divide two dual numbers.
  \<< DC\->R ROT DC\->R
      4. PICK * OVER 4. ROLL * -
      SWAP PICK3 / UNROT SWAP SQ / R\->C
  \>>
  DINV @ Inverse (reciprocal) of a dual number.
  \<< DC\->R OVER INV UNROT
      NEG SWAP SQ / R\->C
  \>>
  DSQ @ Square of a dual number.
  \<< DC\->R OVER SQ UNROT
      * 2. * R\->C
  \>>
  DSQRT @ Square root of a dual number.
  \<< DC\->R SWAP \v/ SWAP OVER 2. * / R\->C
  \>>
  DLN @ Natural logarithm of a dual number.
  \<< DC\->R OVER / SWAP LN SWAP R\->C
  \>>
  DEXP @ Exponentiation (natural antilog) of a dual number.
  \<< DC\->R SWAP EXP SWAP OVER * R\->C
  \>>
  DIPWR @ Dual number to integer power.
  \<< SWAP DC\->R OVER 4. PICK ^
      4. ROLLD SWAP PICK3 1 - ^ * * R\->C
  \>>
  DDPWR @ Dual number to dual number power.
  \<< DC\->R ROT DC\->R OVER 5. PICK ^ 5. ROLLD
      OVER 5. ROLL SWAP OVER 1. - ^ * *
      SWAP LN ROT * PICK3 * + R\->C
  \>>
  DNROOT @ Nth root of dual number for integer n.
  \<< SWAP DC\->R OVER 4. PICK XROOT
      4. ROLLD SWAP ABS PICK3 INV 1. - ^
      * SWAP / R\->C
  \>>
  DC\->R @ Safe C-\>R for either real or dual number.
  \<< DUP TYPE 1. SAME
      :: C\->R 0. IFTE
  \>>
END

Next we have a directory of trig functions including hyperbolic functions, arg and rectangular/polar conversions. They require the command DC→R from the directory above, so that command should be added to this directory if it is not in your path.

Code:


@ Dual number trig functions. All require DC\->R.
DIR
  DSIN
  \<< DC\->R OVER COS *
      SWAP SIN SWAP R\->C
  \>>
  DASIN
  \<< DC\->R SWAP ASIN
      SWAP OVER COS / R\->C
  \>>
  DCOS
  \<< DC\->R OVER SIN NEG *
      SWAP COS SWAP R\->C
  \>>
  DACOS
  \<< DC\->R SWAP ACOS
      SWAP OVER SIN / R\->C
  \>>
  DTAN
  \<< DC\->R SWAP TAN
      SWAP OVER SQ 1 + * R\->C
  \>>
  DATAN
  \<< DC\->R OVER ATAN
      ROT ROT SWAP SQ 1 + / R\->C
  \>>
  DSINH
  \<< DC\->R OVER SINH
      ROT ROT SWAP COSH * R\->C
  \>>
  DASNH
  \<< DC\->R SWAP ASINH
      SWAP OVER COSH / R\->C
  \>>
  DCOSH
  \<< DC\->R OVER COSH
      ROT ROT SWAP SINH * R\->C
  \>>
  DACSH
  \<< DC\->R SWAP ACOSH
      SWAP OVER SINH / R\->C
  \>>
  DTANH
  \<< DC\->R SWAP TANH
      SWAP OVER SQ 1 SWAP - * R\->C
  \>>
  DATNH
  \<< DC\->R OVER ATANH
      ROT ROT SWAP SQ 1 SWAP - / R\->C
  \>>
  DARG @ Argument of dual number.
  \<< DC\->R SWAP /
  \>>
  DR\->P @ Rectangular to polar conversion for dual numbers.
  \<< DC\->R OVER / R\->C
  \>>
  DP\->R @ Polar to Rectangular conversion for dual numbers.
  \<< DC\->R OVER SWAP * R\->C
  \>>
END

Finally, the same directory optimized for the 49/50.

Code:


@ Dual number trig functions. All require DC\->R.
DIR
  DSIN
  \<< DC\->R OVER COS *
      SWAP SIN SWAP R\->C
  \>>
  DASIN
  \<< DC\->R SWAP ASIN
      SWAP OVER COS / R\->C
  \>>
  DCOS
  \<< DC\->R OVER SIN NEG *
      SWAP COS SWAP R\->C
  \>>
  DACOS
  \<< DC\->R SWAP ACOS
      SWAP OVER SIN / R\->C
  \>>
  DTAN
  \<< DC\->R SWAP TAN
      SWAP OVER SQ 1. + * R\->C
  \>>
  DATAN
  \<< DC\->R OVER ATAN
      UNROT SWAP SQ 1. + / R\->C
  \>>
  DSINH
  \<< DC\->R OVER SINH
      UNROT SWAP COSH * R\->C
  \>>
  DASNH
  \<< DC\->R SWAP ASINH
      SWAP OVER COSH / R\->C
  \>>
  DCOSH
  \<< DC\->R OVER COSH
      UNROT SWAP SINH * R\->C
  \>>
  DACSH
  \<< DC\->R SWAP ACOSH
      SWAP OVER SINH / R\->C
  \>>
  DTANH
  \<< DC\->R SWAP TANH
      SWAP OVER SQ 1. SWAP - * R\->C
  \>>
  DATNH
  \<< DC\->R OVER ATANH
      UNROT SWAP SQ 1. SWAP - / R\->C
  \>>
  DARG @ Argument of dual number.
  \<< DC\->R SWAP /
  \>>
  DR\->P @ Rectangular to polar conversion for dual numbers.
  \<< DC\->R OVER / R\->C
  \>>
  DP\->R @ Polar to Rectangular conversion for dual numbers.
  \<< DC\->R OVER SWAP * R\->C
  \>>
END

[Gil]

Very nice and good examples of how to play with User RPN.

And welcome, as new HP50G programs seem not to have (any more?) the real favour of the HP community.

Regards,
Gil

[John Keith]

(02-20-2024 10:17 AM)Gil Wrote:  Very nice and good examples of how to play with User RPN.

Actually they should be rewritten in SysRPL both for improved accuracy using extended reals, and for faster and easier type checking so as to allow exact integer and symbolic versions on the 49 and 50. Unfortunately my SysRPL skills are remedial at best and that is unlikely to happen unless a more competent person wants to volunteer.

[Gil]

As it is now it's most understandable for me and suits me better that way.

[Thomas Klemm]

It took me a while to figure out how to transfer these programs to the iHP48 app.
For this I used the programs IN and OUT that I found in BruceH's post: #2

I copied them, double tapped the screen and used: EDIT → Paste into iHP Stack
Then I manually replaced the trigrams in the string of IN and used OBJ→ to create the program.
This allowed me to use IN to transform OUT.

But you can just as well use these two directly.

IN:

Code:

« →STR 3 TRANSIO RCLF SIZE 3 > # 193357d # 196971d IFTE SYSEVAL + STR→ »

OUT:

Code:

« STD 64 STWS →STR 3 TRANSIO RCLF SIZE 3 > # 193359d # 196297d IFTE SYSEVAL »

Once transferred I moved all programs into a single directory DUAL and created both a CST menu and key assignments in KEYS:

Code:

DIR
  f
    \<< 2 OVER * 3 + OVER DMUL 5 + DMUL 7 +
    \>>
  NEWTON
    \<< 1 R\->C \-> x
      \<< x f C\->R / x C\->R DROP SWAP -
      \>>
    \>>
  CST { { "SINH" { DSINH DASNH } } { "COSH" { DCOSH DACSH } } { "TANH" { DTANH DATNH } } { "R\->P" DR\->P } { "P\->R" DP\->R } { "R\->C" { R\->C DC\->R } } }
  KEYS { S DSIN 41.1 DASIN 41.2 DCOS 42.1 DACOS 42.2 DTAN 43.1 DATAN 43.2 DSQRT 44.1 DSQ 44.2 DNROOT 44.3 DDPWR 45.1 DIPWR 45.2 DARG 45.3 DINV 46.1 DEXP 46.2 DLN 46.3 DDIV 65.1 DMUL 75.1 }
  DSIN
    \<< DC\->R OVER COS * SWAP SIN SWAP R\->C
    \>>
  DASIN
    \<< DC\->R SWAP ASIN SWAP OVER COS / R\->C
    \>>
  DCOS
    \<< DC\->R OVER SIN NEG * SWAP COS SWAP R\->C
    \>>
  DACOS
    \<< DC\->R SWAP ACOS SWAP OVER SIN / R\->C
    \>>
  DTAN
    \<< DC\->R SWAP TAN SWAP OVER SQ 1 + * R\->C
    \>>
  DATAN
    \<< DC\->R OVER ATAN ROT ROT SWAP SQ 1 + / R\->C
    \>>
  DSINH
    \<< DC\->R OVER SINH ROT ROT SWAP COSH * R\->C
    \>>
  DASNH
    \<< DC\->R SWAP ASINH SWAP OVER COSH / R\->C
    \>>
  DCOSH
    \<< DC\->R OVER COSH ROT ROT SWAP SINH * R\->C
    \>>
  DACSH
    \<< DC\->R SWAP ACOSH SWAP OVER SINH / R\->C
    \>>
  DTANH
    \<< DC\->R SWAP TANH SWAP OVER SQ 1 SWAP - * R\->C
    \>>
  DATNH
    \<< DC\->R OVER ATANH ROT ROT SWAP SQ 1 SWAP - / R\->C
    \>>
  DARG
    \<< DC\->R SWAP /
    \>>
  DR\->P
    \<< DC\->R OVER / R\->C
    \>>
  DP\->R
    \<< DC\->R OVER SWAP * R\->C
    \>>
  DINV
    \<< DC\->R OVER INV ROT ROT NEG SWAP SQ / R\->C
    \>>
  DMUL
    \<< DC\->R ROT DC\->R OVER 5 PICK * 5 ROLLD ROT ROT * ROT ROT * + R\->C
    \>>
  DDIV
    \<< DC\->R ROT DC\->R 4 PICK * OVER 4 ROLL * - SWAP 3 PICK / ROT ROT SWAP SQ / R\->C
    \>>
  DSQ
    \<< DC\->R OVER SQ ROT ROT * 2 * R\->C
    \>>
  DSQRT
    \<< DC\->R SWAP \v/ SWAP OVER 2 * / R\->C
    \>>
  DLN
    \<< DC\->R OVER / SWAP LN SWAP R\->C
    \>>
  DEXP
    \<< DC\->R SWAP EXP SWAP OVER * R\->C
    \>>
  DIPWR
    \<< SWAP DC\->R OVER 4 PICK ^ 4 ROLLD SWAP 3 PICK 1 - ^ * * R\->C
    \>>
  DDPWR
    \<< DC\->R ROT DC\->R OVER 5 PICK ^ 5 ROLLD OVER 5 ROLL SWAP OVER 1 - ^ * * SWAP LN ROT * 3 PICK * + R\->C
    \>>
  DNROOT
    \<< SWAP DC\->R OVER 4 PICK INV ^ 4 ROLLD SWAP ABS 3 PICK INV 1 - ^ * SWAP / R\->C
    \>>
  DC\->R
    \<< DUP TYPE 1 SAME
      \<< C\->R
      \>> 0 IFTE
    \>>
END

Once the directory object has been copied over and translated by IN just save it in a variable say DUAL.
This will create the directory with all the files.

The contents of KEYS can then be stored using STOKEYS.

Then hit USER twice to activate the assignments permanently.
Use CST to activate the customer menu.

Examples:

Arithmetic

\(
(3 + 4\varepsilon)(5 + 6\varepsilon)
\)

3 4 R→C
5 6 R→C
*

(15,38)

Evaluation of a function and its derivative

\(
\begin{align}
f(x) = \frac{1}{\sqrt{3 + \sin(x)}}
\end{align}
\)

Evaluate \(f(2)\) and \(f'(2)\).

2 1 R→C
SIN
3 +
√
INV

(0.505767179164,2.69195956021E-2)

Polynomial

Write a program to evaluate the following polynomial:

\(
f(x) = 2x^3 + 3x^2 + 5x + 7
\)

Code:

\<< 
  2 OVER *
  3 + OVER DMUL
  5 + DMUL
  7 +
\>>

Hint: Please note the use of the ordinary multiplication × with a constant.

Finding a root with Newton's algorithm

\(
\begin{align}
x \mapsto x - \frac{f(x)}{f'(x)}
\end{align}
\)

The program is straight forward since both the function and its derivative can be calculated with a single call to f:

Code:

\<<
  1 R\->C \-> x
  \<<
    x f C\->R /
    x C\->R
    DROP SWAP -
  \>>
\>>

Now we can use it to find the root of the previous polynomial:

-1 NEWTON NEWTON NEWTON …

-1.6
-1.4594795539
-1.44565202879
-1.44552846845
-1.44552845868

-1.44552845868

[Gil]

The following precious program INOUT was once given to me by one of the HP Museum member:

Code:

\<< RCWS RCLF \-> ws f
  \<< 3 TRANSIO DUP TYPE 2
    IF ==
    THEN \->STR f SIZE 3 > # 2F34Dh # 3016Bh IFTE SYSEVAL + STR\->
    ELSE STD 64 STWS \->STR f SIZE 3 > # 2F34Fh # 2FEC9h IFTE SYSEVAL
    END ws STWS f STOF
  \>>
\>>

[Thomas Klemm]

(02-21-2024 04:49 AM)Gil Wrote:  The following precious program INOUT was once given to me by one of the HP Museum member:

It's from post #4 by user grsbanks in the same thread I linked earlier.

[DavidM]

(02-20-2024 10:17 AM)Gil Wrote:  And welcome, as new HP50G programs seem not to have (any more?) the real favour of the HP community.

I don't think this is accurate. It seems to me that there are a handful of vocal members here at the hpmuseum who seem to enjoy disparaging the RPL calculators, and I'd even go so far as to say that there has been some small amount of "chilling effect" regarding participation in RPL-related threads as a result. But I don't believe that should be interpreted as a consensus of the community's appreciation for those systems.

Case in point: Eddie Shore recently posted a poll regarding preferred systems for "quick number crunching", and the RPL systems were well-represented in the results.

Personally, I find threads like this one intriguing. I just lack the mathematical background to make practical use of the concepts. I'd enjoy participating in an effort to implement these commands in SysRPL, though (time permitting). Certain characteristics of the 49g-50g SysRPL environment are a good fit for this kind of application. In particular, I'd like to better understand the implications of function overloading for exact integers here vs. extended reals.

[Gil]

I meant many RPN programs for calculators, yes, but for other calculators than the HP50G.

Just an impression.

Regards

[Thomas Klemm]

(02-21-2024 04:34 PM)DavidM Wrote:  Personally, I find threads like this one intriguing. I just lack the mathematical background to make practical use of the concepts.

Consider two resistors with a certain tolerance:

\(
\begin{align}
R_1 &= 3 \pm 4 \varepsilon \\
R_1 &= 5 \pm 6 \varepsilon \\
\end{align}
\)

What is the total resistance and its tolerance if both resistors are in parallel?

3 4 R→C 1/X
5 6 R→C 1/X
+ 1/X

(1.875,2.40625)


Compare this to the calculation if \(\varepsilon = 1e-5\):

3.00004 1/X
5.00006 1/X
+ 1/X

1.8750240625


With each number we include the linearization of its perturbation.

[John Keith]

Next up, dual number elliptic integrals... and a mystery.

The programs DNEK, DNEE, DNEΠ and DNEN compute respectively the dual number versions of K(k), E(k), Π(n, k) and q(k), the elliptic nome. They use the programs from this thread on elliptic integrals. I am also including translations of the programs DAGM and DELK from the Double ROM manual.

Á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. Perhaps someone who is an HP-41 user familiar with the Double ROM can check my results and explain the different results.

I am listing the program as a directory including all dual number and elliptic integral programs needed to run them. The directory is thus self-contained and can be used with any RPL calculator, although the names of programs for Π(n, k) will have to be changed for the HP-28 due to its character set not having Π.

Code:


DIR
  DNEK
  \<< DC\->R OVER Kk
      ROT ROT SWAP dKk * R\->C
  \>>
  DNEE
  \<< DC\->R OVER EKk DROP
      ROT ROT SWAP dEk * R\->C
  \>>
  DNE\PI
  \<< DC\->R ROT ROT DUP2 d\PIdk
      ROT ROT \PInk ROT ROT * R\->C
  \>>
  DNEN
  \<< DC\->R OVER ENOME
      ROT ROT SWAP dNOME * R\->C
  \>>
  DAGM
  \<<
    IF DUP RE
    THEN
      DO DUP ROT ROT DUP2 + 2 DDIV
        ROT ROT DMUL DSQRT
      UNTIL ROT OVER ==
      END DROP
    END
  \>>
  DELK
  \<< \pi \->NUM 1 ROT DSQ - DSQRT
      1 SWAP DAGM DUP + DDIV
  \>>
  dKk
  \<< 1 OVER SQ - OVER EKk
      3 PICK * - ROT ROT * /
  \>>
  dEk
  \<< DUP EKk - SWAP /
  \>>
  dNOME
  \<< DUP ENOME \pi SQ \->NUM *
    SWAP 1 OVER SQ - OVER Kk SQ
    ROT 2 * * * /
  \>>
  d\PIdk
  \<<
    IF OVER ABS
    THEN DUP2 \PInk ROT ROT DUP EKk DROP
      OVER SQ 1 - /
      4 ROLL + OVER *
      ROT ROT SQ - /
    ELSE NIP dKk
    END
  \>>
  Kk
  \<< DUP ABS 1
    IF \=/
    THEN 1 1 ROT SQ - \v/
      DO DUP ROT ROT DUP2 + 2 /
        ROT ROT * \v/
      UNTIL ROT OVER ==
      END DROP 2 * \pi \->NUM SWAP /
    ELSE DROP MAXR \->NUM
    END
  \>>
  EKk
  \<< DUP ABS 1
    IF \=/
    THEN 1 SWAP SQ DUP2 - \v/
      ROT 3 PICK 2 / 4 ROLLD .5 ROT ROT
      WHILE ROT 2 * ROT ROT
        DUP ROT ROT DUP2 * \v/
        ROT ROT + 2 / ROT OVER \=/
      REPEAT 4 ROLL OVER 4 * / SQ
        4 PICK OVER * 6 ROLL +
        5 ROLLD 4 ROLLD
      END 5 ROLLD DROP2 DROP
      \pi \->NUM SWAP 1 SWAP - OVER *
      3 PICK 2 * / ROT 2 * ROT SWAP /
    ELSE DROP 1 MAXR \->NUM
    END
  \>>
  PInk
  \<<
    IF OVER ABS
    THEN
      IF OVER 1 == OVER ABS 1 == OR
      THEN DROP2 MAXR \->NUM
      ELSE
        IF OVER DUP RE 1 > SWAP IM NOT AND
        THEN DROP2 # 203h DOERR
        ELSE SWAP 1 1 \-> n s q
          \<< 1 1 ROT SQ - \v/ OVER n - ROT ROT
            DO 3 DUPN * DUP2 + ROT ROT - OVER /
              q * 2 / s 6 ROLLD 6 PICK
              OVER + 's' STO 'q' STO
              4 ROLL \v/ 2 * / SQ ROT ROT
              DUP ROT ROT DUP2 + 2 / ROT ROT * \v/
            UNTIL ROT OVER == 5 ROLL s == AND
            END ROT ROT DROP2
            s n 1 OVER - / * 2 +
            \pi \->NUM * SWAP 4 * /
          \>>
        END
      END
    ELSE NIP Kk
    END
  \>>
  ENOME
  \<< 1 OVER SQ - \v/ Kk
      \pi NEG \->NUM * SWAP Kk / EXP
  \>>
  DC\->R
  \<< DUP TYPE 1 SAME
    \<< C\->R
    \>> 0 IFTE
  \>>
  DMUL
  \<< DC\->R ROT DC\->R
      OVER 5 PICK * 5 ROLLD
      ROT ROT * ROT ROT * + R\->C
  \>>
  DDIV
  \<< DC\->R ROT DC\->R
      4 PICK * OVER 4 ROLL * -
      SWAP 3 PICK / ROT ROT SWAP SQ / R\->C
  \>>
  DSQ
  \<< DC\->R OVER SQ
      ROT ROT * 2 * R\->C
  \>>
  DSQRT
  \<< DC\->R SWAP \v/
      SWAP OVER 2 * / R\->C
  \>>
END

And a listing optimized for the HP 49 and 50.

Code:


DIR
  DNEK
  \<< DC\->R OVER Kk
      UNROT SWAP dKk * R\->C
  \>>
  DNEE
  \<< DC\->R OVER EKk DROP
      UNROT SWAP dEk * R\->C
  \>>
  DNE\PI
  \<< DC\->R UNROT DUP2 d\PIdk
      UNROT \PInk UNROT * R\->C
  \>>
  DNEN
  \<< DC\->R OVER ENOME
      UNROT SWAP dNOME * R\->C
  \>>
  DAGM
  \<<
    IF DUP RE
    THEN
      DO DUP UNROT DUP2 + 2. DDIV
        UNROT DMUL DSQRT
      UNTIL ROT OVER ==
      END DROP
    END
  \>>
  DELK
  \<< \pi \->NUM 1. ROT DSQ - DSQRT
      1. SWAP DAGM DUP + DDIV
  \>>
  dKk
  \<< 1. OVER SQ - OVER EKk
    PICK3 * - UNROT * /
  \>>
  dEk
  \<< DUP EKk - SWAP /
  \>>
  dNOME
  \<< DUP ENOME \pi SQ \->NUM *
    SWAP 1. OVER SQ - OVER Kk SQ
    ROT 2. * * * /
  \>>
  d\PIdk
  \<<
    IF OVER ABS
    THEN DUP2 \PInk UNROT DUP EKk DROP
      OVER SQ 1. - /
      4. ROLL + OVER *
      UNROT SQ - /
    ELSE NIP dKk
    END
  \>>
  Kk
  \<< DUP ABS 1.
    IF \=/
    THEN 1. 1. ROT SQ - \v/
      DO DUP UNROT DUP2 + 2. / UNROT * \v/
      UNTIL ROT OVER ==
      END DROP 2. * \pi \->NUM SWAP /
    ELSE DROP MAXR \->NUM
    END
  \>>
  EKk
  \<< DUP ABS 1.
    IF \=/
    THEN 1. SWAP SQ DUP2 - \v/
      ROT PICK3 2. / 4. ROLLD .5 UNROT
      WHILE ROT 2. * UNROT
        DUP UNROT DUP2 * \v/
        UNROT + 2. /
        ROT OVER \=/
      REPEAT 4. ROLL OVER 4. * / SQ
        4. PICK OVER * 6. ROLL +
        5. ROLLD 4. ROLLD
      END 5. ROLLD DROP2 DROP
      \pi \->NUM SWAP 1. SWAP - OVER *
      PICK3 2. * / ROT 2. * ROT SWAP /
    ELSE DROP 1. MAXR \->NUM
    END
  \>>
  \PInk
  \<< IF OVER ABS
      THEN
        IF OVER 1. == OVER ABS 1. == OR
        THEN DROP2 MAXR \->NUM
        ELSE
          IF OVER DUP RE 1. > SWAP IM NOT AND
          THEN DROP2 #203h DOERR
          ELSE SWAP 1. 1. \-> n s q
            \<< 1. 1. ROT SQ - \v/ OVER n - UNROT
              DO PICK3 PICK3 PICK3 *
                DUP2 + UNROT - OVER /
                q * 2. / s 6. ROLLD 6. PICK
                OVER + 's' STO 'q' STO
                4. ROLL \v/ 2. * / SQ UNROT
                DUP UNROT DUP2 + 2. / UNROT * \v/
              UNTIL ROT OVER == 5. ROLL s == AND
              END UNROT DROP2
              s n 1. OVER - / * 2. +
              \pi \->NUM * SWAP 4. * /
            \>>
          END
        END
      ELSE NIP Kk
      END
  \>>
  ENOME
  \<< 1. OVER SQ - \v/ Kk
      \pi NEG \->NUM * SWAP Kk / EXP
  \>>
  DC\->R
  \<< DUP TYPE 1. SAME :: C\->R 0. IFTE
  \>>
  DMUL
  \<< DC\->R ROT DC\->R OVER 5. PICK * 5. ROLLD
      UNROT * UNROT * + R\->C
  \>>
  DDIV
  \<< DC\->R ROT DC\->R 4. PICK * OVER 4. ROLL * -
      SWAP PICK3 / UNROT SWAP SQ / R\->C
  \>>
  DSQ
  \<< DC\->R OVER SQ UNROT * 2 * R\->C
  \>>
  DSQRT
  \<< DC\->R SWAP \v/ SWAP OVER 2 * / R\->C
  \>>
END

[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

[Albert Chan]

It may be simpler to use elliptic parameter m = k^2
Example, here is prove for Legendre's relation using m

https://en.wikipedia.org/wiki/Elliptic_integral Wrote:Complete elliptic integral of the second kind, Derivative and differential equation

With f' = df/dm, m^* = 1-m, above differential equations simplified to:

E' = (E-K) / (2 m)
E = -4 m^* (m E')' = -2 m^* (E-K)'

-E / (2 m^*) = E' - K'
K' = (E-K) / (2 m) + E / (2 m^*)
2 m m^* K' = m^* (E-K) + m E = E - m^* K = (E-K) + m K
(E-K) = 2 m m^* K' - m K = (2 m) (m^* K' - K/2)

E' = m^* K' - K/2
E = K + 2 m E'

Redo previous post example using m = k^2

lua> D = require'dual'.D
lua> m = D.new(0.5^2, 1)
lua> a, b = 1, (1-m):sqrt()
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( m=0.25+ε )
1.685750354812596        0.5417318486132805

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

Here, we get same numbers by chance, because dm/dk = 2k = 2*0.5 = 1

[John Keith]

Thanks for confirming my calculations. I was a bit hesitant to question the Double ROM manual and its esteemed author.

I considered using m instead of k, and I may yet reconsider as using m does make some things easier.