Extended Complex Number Calculations in SysRPL


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?



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!



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.


Possibly Related Threads…
Thread Author Replies Views Last Post
  HP Prime: complex numbers in CAS. Alberto Candel 1 1,895 12-06-2013, 02:36 PM
Last Post: parisse
  [HP Prime] Plots containing complex numbers bug? Chris Pem10 7 3,622 12-05-2013, 07:40 AM
Last Post: cyrille de Brébisson
  Complex Number Entry on Prime Jeff O. 19 5,121 11-16-2013, 12:34 PM
Last Post: Jeff O.
  HP Prime complex results Javier Goizueta 0 976 10-06-2013, 12:59 PM
Last Post: Javier Goizueta
  HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 4,212 09-30-2013, 03:44 AM
Last Post: From Hong Kong
  [HP-Prime xcas] operations with complex numbers + BUGs + Request CompSystems 9 3,431 09-08-2013, 10:40 PM
Last Post: CompSystems
  ILPer with "auto-extended addressing" Christoph Giesselink 0 978 07-23-2013, 03:28 PM
Last Post: Christoph Giesselink
  Elliptic integrals of 1st and 2nd kind calculated by complex agm Gjermund Skailand 3 1,455 06-29-2013, 03:39 PM
Last Post: Gjermund Skailand
  HP-11C and complex numbers Antlab 5 2,194 06-28-2013, 08:59 AM
Last Post: Antlab
  71B Complex User Defined Functions in Calc Mode? Michael Burr 0 927 03-18-2013, 09:38 PM
Last Post: Michael Burr

Forum Jump: