(28/48/50) Dual Number Functions
02-19-2024, 08:37 PM
Post: #1
 John Keith Senior Member Posts: 1,027 Joined: Dec 2013
(28/48/50) Dual Number Functions
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
02-20-2024, 10:17 AM (This post was last modified: 02-20-2024 10:27 AM by Gil.)
Post: #2
 Gil Senior Member Posts: 586 Joined: Oct 2019
RE: (28/48/50) Dual Number Functions
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
02-20-2024, 09:23 PM
Post: #3
 John Keith Senior Member Posts: 1,027 Joined: Dec 2013
RE: (28/48/50) Dual Number Functions
(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.
02-21-2024, 12:57 AM
Post: #4
 Gil Senior Member Posts: 586 Joined: Oct 2019
RE: (28/48/50) Dual Number Functions
As it is now it's most understandable for me and suits me better that way.
02-21-2024, 04:34 AM
Post: #5
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: (28/48/50) Dual Number Functions
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
02-21-2024, 04:49 AM
Post: #6
 Gil Senior Member Posts: 586 Joined: Oct 2019
RE: (28/48/50) Dual Number Functions
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   \>> \>>
02-21-2024, 05:56 AM
Post: #7
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: (28/48/50) Dual Number Functions
(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.
02-21-2024, 04:34 PM
Post: #8
 DavidM Senior Member Posts: 982 Joined: Dec 2013
RE: (28/48/50) Dual Number Functions
(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.
02-21-2024, 04:47 PM
Post: #9
 Gil Senior Member Posts: 586 Joined: Oct 2019
RE: (28/48/50) Dual Number Functions
I meant many RPN programs for calculators, yes, but for other calculators than the HP50G.

Just an impression.

Regards
02-22-2024, 07:13 AM
Post: #10
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: (28/48/50) Dual Number Functions
(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.
03-20-2024, 07:14 PM (This post was last modified: 03-20-2024 07:16 PM by John Keith.)
Post: #11
 John Keith Senior Member Posts: 1,027 Joined: Dec 2013
RE: (28/48/50) Dual Number Functions
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
03-21-2024, 05:13 PM (This post was last modified: 03-22-2024 07:49 PM by Albert Chan.)
Post: #12
 Albert Chan Senior Member Posts: 2,551 Joined: Jul 2018
RE: (28/48/50) Dual Number Functions
(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

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
03-22-2024, 04:23 PM (This post was last modified: 03-22-2024 07:50 PM by Albert Chan.)
Post: #13
 Albert Chan Senior Member Posts: 2,551 Joined: Jul 2018
RE: (28/48/50) Dual Number Functions
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
03-22-2024, 07:04 PM
Post: #14
 John Keith Senior Member Posts: 1,027 Joined: Dec 2013
RE: (28/48/50) Dual Number Functions
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.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)