[wp34s] Yet another complex solver

[emece67]

Hi all. Having some trouble using the complex solver in the wp34s library, I've developed a new complex solver to be (hopefully) included in the wp34s library.

In my case, the complex solver already in the library does not converge many times. I suppose this is some kind of glitch on its wp34s porting (the original program was for the 35s), as the article by its author (Valentín Albillo) describing the program states 5 example applications some of them are not solved by the wp34s incarnation.

Anyway, I wrote a new complex solver for the wp34s. In my tests it behaves quite correctly, is fast and accurate. It is based upon the Ostrowski-Traub method, a 4th order convergent algorithm (Newton is 2nd order, Laguerre 3rd) that only needs 3 function evaluations on each iteration.

In fact it is a heavy rewriting of another complex solver I wrote years ago to be used on my first HP calculator, an hp28C. I used such solver when studying EE to solve many waveguide problems.

I paste here the header of the programs listing, which describes its usage and some details of the algorithm. I've called it CSV and occupies 161 program steps. Can be used in both DBLOFF & DBLON modes. Sometimes it is even faster than the built in (real) solver when solving for real roots of real valued functions.

Code:


/*
* CSV, an Ostrowski-Traub based complex solver
*
*   v1.4 (20151222) by M. César Rodríguez GGPL
*
* Usage:
*   Input:
*     XY: contains an initial root estimation. X is the real part, Y the
*         imaginary one
*     alpha: contains the name of the program to be solved. This program
*         is called with the stack (an SSIZE8 one) filled with a complex
*         estimation of the root. The program must return a complex
*         number in XY (X real, Y imaginary)
*
*   Output:
*     XY: the root found or the value of the last estimation for the
*         root. X real, Y imaginary
*     ZT: the value of the program to be solved evaluated at the
*         returned complex value in XY. Z real, T imaginary
*     LI: the input XY value
*
*     If the user interrupts the program ([<-] key), or the program to
*       be solved evaluates to infinity or NaN, or the solver reaches
*       its maximum iteration number (800 iterations), it ends with the
*       message "No root Found"
*
*   The user can interrupt the program pressing the [<-] key. If the
*     user interrupts the program with [EXIT] or [R/S], the values in
*     the stack, ABCD, L & I registers and flag D are not defined. Also,
*     the size of the stack will be 8
*
*   Register usage:
*     ABCD are cleared (no matter if the user works in SSIZE4 or SSIZE8)
*     Stack size is forced to 8 levels, but it is restored to the value
*       prior to calling the solver upon program end
*     Flag D is set (so function to be solved returns infinite or NaN
*       instead of reporting errors). This flag is restored to the value
*        prior to calling the solver upon program end
*
* Algorithm:
*   The Ostrowski-Traub algorithm has 4th order convergence into
*     single roots. Requires 3 function evaluations on each iteration.
*     See, for example, [1]
*
*   For each root estimate z_i, a new estimate is computed with:
*       w_i = z_i - f(z_i)/f'(z_i)    (this is a Newton iteration)
*
*     where here the derivative f'(z_i) is estimated with:
*       f'(z_i) ~ (f(z_i + e_i) - f(z_i))/e_i
*
*   Thus, the expression used here to get the w_i estimate is:
*       w_i = z_i + e_i/(1 - f(z_i + e_i)/f(z_i))
*
*   Then, such estimate is refined with:
*       z_{i+1} = w_i - f(w_i)(z_i - w_i)/(f(z_i) - 2f(w_i))
*
*     which  is implemented here as:
*       z_{i+1} = w_i + (z_i - w_i)/(2 - f(z_i)/f(w_i))
*
*   After each iteration, |z_{i+1} - z_i| is checked to be smaller
*     than max(|z_i|/10^4, 10^-8). If so, a flag is set signaling that
*     the algorithm is "converging". This flag remains set until
*     program end
*
*   After each iteration and if the algorithm is "converging",
*     |z_{i+1} - z_i| is compared to |z_i - z_{i-1}|, being greater the
*     algorithm is supposed to had converged and z_i (not z_{i+1}) is
*     declared a root. This stopping algorithm is, basically, the Ward's
*     criterion [2]
*
*   The program will also exit with a root found (hopefully) if:
*     - evaluating the program to be solved returns 0 + 0j
*     - f(z_i + e_i) = f(z_i)
*     - z_{i+1} = z_i
*
*   The program will exit with error message "No root Found" if:
*     - evaluating the program to be solved returns Infinite or NaN
*     - the algorithm reaches 800 iterations (in DBLOFF mode and with
*         single roots, the algorithm many times needs only 3 iterations
*         to converge)
*     - the user interrupts the program pressing [<-]
*
*   When a root is found, it is "beautified" by:
*     - changing numbers like 1.124999999973587 or 1.125000000032518
*       to 1.125. The code for this is adapted from the PRS polynomial
*       solver for the wp34s by Franz Huber
*     - numbers with an imaginary (real) part much smaller than its
*       real (imaginary) part, are converted to real (imaginary) ones
*       zeroing the smaller part
*     Then, the function to be solved is evaluated at this new,
*     beautified, root and it the value of the function is closer
*     (or equally close) to 0, then the beautified root is returned as
*     the root
*
*   When estimating the value of the derivative f'(z_i) with:
*     (f(z_i + e_i) - f(z_i))/e_i, the increment e_i is:
*     - if z_i is real, then
*         e_i = |z_i - z_{i-1}|/10^3
*     - otherwise
*         e_i = |z_i - z_{i-1}|/10^3·exp(j·ai), with a_i a random
*         angle between 0 and 2·pi
*     Thus, if the function to be solved is real valued for real
*     arguments and the initial root estimate is real, the program
*     will only search for real roots
*
* When called from another program, this solver will behave like a
*   test, executing next program step when a root is found, or
*   skipping it when not found
*
* This solver does also work in DBLON mode
*
* References:
*   [1] Wei-hong, B.I., Qing-biao, W.U., Hong-min R. "Convergence ball
*     and error analysis of Ostrowski-Traub’s method." Appl. Math. J.
*     Chinese Univ. 2010, 25(3): 374-378
*
*   [2] Ward, R.C. "The QR algorithm and Hyman’s method on vector
*     computers." Mathematics of Computation 01/1976; 30(133): 132-132
*
*/

Find attached the full listing. Any bug reporting or suggestion will be welcome.

EDIT: last version moved to the General Software Library here

Regards.

[fhub]

(12-23-2015 02:21 AM)emece67 Wrote:  Find attached the full listing. Any bug reporting or suggestion will be welcome.

WOW, really a great program - and what a detailled description!

When you first wrote about it some time ago, I thought you would talk about a complex polynomial solver, but now it's even a 'full' complex solver -
it certainly deserves to be included in the official WP34s library!

Many thanks,
Franz.

[Marcus von Cube]

(12-23-2015 05:01 PM)fhub Wrote:  it certainly deserves to be included in the official WP34s library!

Chances are high that it will be included not only in source form but also in the _full images if no-one objects.

[walter b]

I concur.

d:-)

[emece67]

(12-23-2015 05:01 PM)fhub Wrote:  When you first wrote about it some time ago, I thought you would talk about a complex polynomial solver, but now it's even a 'full' complex solver

In fact I was talking about a poly solver, but after your superb PRS program, I found not a big motivation to continue such way, so I switched to the general complex solver.

I'm thinking now in adding support for polynomial operations to the wp34s (add, substract, multiply, divide, power, from roots get poly an root finding) to be added to the MATRIX menu. It would be an addition to the Matrix Reloaded patches I uploaded a year ago.

(12-23-2015 06:39 PM)Marcus von Cube Wrote:  Chances are high that it will be included not only in source form but also in the _full images if no-one objects.

If so, please, use the CSV.zip I've uploaded a moment ago. It has a minor enhancement.

[Marcus von Cube]

(12-24-2015 09:59 AM)emece67 Wrote:  If so, please, use the CSV.zip I've uploaded a moment ago. It has a minor enhancement.

It will appear soon on SF. BTW, I added the poly root solver to the compiled library at the same time but had to remove the latter from the calc_ir_full build because flash is full. :-(

I did not touch the "complex mode" branch which I'm not familiar enough with.