HP Forums
Extended Complex Number Calculations in SysRPL - Printable Version

+- HP Forums (https://archived.hpcalc.org/museumforum)
+-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html)
+--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html)
+--- Thread: Extended Complex Number Calculations in SysRPL (/thread-117286.html)



Extended Complex Number Calculations in SysRPL - Les Wright - 06-25-2007

This follows up a back and forth I was having with James Prange in an earlier thread.

I have completed my SysRPL code to compute the sine and cosine integrals by complex continued fraction. The original Modified Lentz algorithm and formulae can be found in Numerical Recipes 2nd ed or later. Internally, the complex numbers are treated as two part entities, where each part is an extended real. The parts can be either the real or imaginary parts of rectangular form, or the modulus and argument of polar form. I have used NULLLAMs and stack manipulation. The complex arithmetic really is normal extended real number arithmetic on the constituent parts--I work in rectangular form for additions and subtractions, and polar form for multiplication, division, roots, and powers. The SysRPL functions %%R>P and %%P>R allow ready switch back and forth without getting bogged down in the trigonometry.

I have managed to keep track of everything and it works fine, but this is one case where the SysRPL code, albeit not yet optimized is perceptibly slower than my UserRPL version. However, the trade of can be construed as worth it--since the internal calculations are done to 15 digits, the 12 digit output is completely accurate provided that the input argument is of suitable size.

This leads me to a practical question, that I believe someone asked here earlier. When working in extended reals in SysRPL, what would you chose as a "machine epsilon" for the purpose of testing convergence? In double precision on my PC, I think machine epsilon is about 2.2e-16. In my program, the convergence test is when a particular multiplicative factor becomes unity and doesn't change in a subsequent iteration. In UserRPL, converging to actual equality works fine, but in extended reals in SysRPL, this expectation drags the computation out unnecessarily for some input. I have by trial and error settled on 2e-15 as my EPS. Any smaller and you might as well just converge to equality unnecessarily. So if variable del is the factor I want to be unity, the loop exits when LAM del %%1 %%- %%ABS %%2E-15 %%< returns True. Requiring LAM del %%1 EQUAL to be true consumes cpu time and offers no added benefit.

Any thoughts on this?

Les




Re: Extended Complex Number Calculations in SysRPL - Les Wright - 06-26-2007

There is potential for significant speed improvement here.

Though Dave has yet to add it to the HP41 Software Library, JM Baillard has shared with me a routine the computes the continued fraction from right to left by brute force for a predetermined number of terms. There is no futzing around with convergence tests, and flopping back and forth between polar and rectangular forms of the complex numbers is kept to a minimum. It is easy enough to determine the minimum argument size that will yield full precision for a given number of terms, and above that the extra terms are computed unnecessarily, but in a fast environment this shouldn't be that noticeable. Certainly it should be less noticeable than the delays in my own code due to convergence tests and frequent changes between polar and rectangular mode. I should try to port that HP41 routine (which is, btw, less than 50 steps in length) and see how I make out. This may be one case where using the Modified Lentz algorithm is like using a shotgun to kill a fly!

Les


Re: Extended Complex Number Calculations in SysRPL (JMBaillard take note!) - Les Wright - 06-29-2007

My port of JMB's code was a success, but it still isn't as fast as I would like. I gives results in around a second to full 12 digit precision for arguments greater than 5 (arguments less than 5 are well treated by a series summation which proceeds much faster). I may post my code here and in comp.sys.hp48 to see if I can get some optimization advice.


Les