Post Reply 
(34S) Real/complex solver
04-16-2017, 05:28 PM (This post was last modified: 06-15-2017 01:21 PM by Gene.)
Post: #1
(34S) Real/complex solver
This is an upgrade to my previous Ostrowsky-Traub complex solver for the wp-34s. Now it can also be used to, directly, solve real equations.

Call is 'CSV' to solve complex equations, call it 'RSV' and it will solve real ones.

I've also modified some internal parameters so it is now slightly more precise. The root "beautifier" has also been upgraded, now it looks, in the final root estimation, for strings of four or more 0s or 9s and rounds the estimation accordingly. The beautified root is checked (evaluating the equation at it) to be a better estimation than the bare one.

As a bonus, each estimation (real or complex) is displayed on the screen, so you can see the algorithm progress (as the built-in integration program does).

There are two versions of this program:
  • v1.6R, Uses local labels and is intended to be entered by hand
  • v1.6L, uses long labels and is intended to be assembled

Comparing it with the (real) built-in solver, is needs sometimes less iterations to get a root, and more in some other times, not a clear difference here, but as this program works internally in complex mode, it uses to be slower than the built-in solver. On the other side, it seems to be more robust than the built-in (fails to get a root in less cases), shows all estimations in the display (you can see how the algorithm progresses) and only needs an initial estimation of the root (not a interval, as the built-in needs. This, for me, is a benefit).

Here is the v1.6R:
Code:

/*
* RCSV, an Ostrowski-Traub based real/complex solver
*
*   v1.6R (20170414) by M. César Rodríguez GGPL
*
* Usage:
*   When used as a real solver (called as 'RSV'):
*     Input:
*       X:      contains an initial root estimation
*       alpha:  contains the name of the program to be solved. This program
*               is called with the stack (an SSIZE8 one) filled with an
*               estimation of the root. The program must return a number in X
*     Output:
*       X:      the found root or the value of the last estimation for the root
*       Y:      the value of the program to be solved evaluated at the returned
*               value in X
*       ZT:     lost (zeroed)
*       L:      the input X value
*       IJK:    kept intact
*       ABCD:   lost (zeroed)
*
*   When used as a complex solver (called as 'CSV'):
*     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
*       JK:     kept intact
*       ABCD:   lost (set to zero)
*
*     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 (400 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. The same will occur if the
*     program to be solved throws an error (e.g. trying to access a non
*     existent register)
*
*   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 for all internal calculations in order not to throw
*       errors when a calculation returns Infinite or NaN, but the program
*       to be solved is called with the flag D the user set. The D flag
*       is also restored at 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 approximated 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^-8, 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 (0 in real mode)
*     - 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 400 iterations
*     - the user interrupts the program pressing [<-]
*
*   When a root is found, it is "beautified" by:
*     - changing numbers like 1.124999973587 or 1.125000032518
*       to 1.125. The code looks for chains of four or more 0s or 9s and
*       rounds the number accordingly
*     - in complex mode, 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, in complex mode, 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
*
* 230 steps, 15 local registers
*
*/
        LBL'RSV'            // Real SolVer
            LocR 15
            SF .03          // flag the real mode
            GTO 07

        LBL'CSV'            // Complex SolVer
            LocR 15         // Local registers and flags:
                            //  .00       mode bits
                            //  .01       |d_i| (that's |z_i - z_{i-1}|)
                            //  .02-03    z_i
                            //  .04-05    e_i
                            //  .06-07    f(z_i)
                            //  .08-09    z_0
                            //  .10-11    arg, argument for program to
                            //            solve, also beautified zi
                            //  .12-13    w_i
                            //  .14       iteration counter
                            //  F.00      no convergence, or f evaluates
                            //            to special, or user interrupt
                            //            ([<-])
                            //  F.01      algorithm is converging
                            //  F.02      backup of D flag
                            //  F.03      true in real mode

          LBL 07            // common entry
            STOM .00        // STOM mode. Mode is restored at exit
            SSIZE8
            FS? D           // backup flag D
              SF .02
            SF D            // specials are not errors here
            FS? .03         // in real mode, I must be maintained, so
              y<> I         // it is saved as the imaginary part of Z_0
            cSTO .08        // cSTO z0. Initial guess
            FS? .03         // in real mode, the imaginary part of the
              # 000         // estimates should be 0
            FS? .03
              x<> Y
            cSTO .02        // cSTO zi. Estimate
            cABS            // |z0|
            x=0?            // |z0| = 0?
              # 001
            SDR 003         // |d0| = |z0|*10-3, or 10^-3 if |z0| = 0
            STO .01         // later it is divided again by 10^3. STO |di|
            # 040
            SDL 001         // 400, maximum iter #
            STO .14         // STO counter
            cRCL .02        // cRCL zi

            // first pass, Newton's estimate (zi on stack)
          LBL 14            // main loop
            cSTO .02        // cSTO zi
            TOP?            // show progress
              XEQ 91
            XEQ 84          // eval function, check for special or root
              GTO 49        // if result is special or root, exit
            cSTO .06        // cSTO f(zi)

            // compute e_i to estimate derivative
            RCL .03         // RCL Im(zi)
            # 000
            x=? Y           // zi is real?
              GTO 21        // yes, e_i must also be real
            ACOS            // no. pi/2 (or 90º or 100 grad)
            # 004
            *
            RAN#
            *               // ai, random angle between 0 and 2pi
          LBL 21
            RCL .01         // RCL |di|. ai or 0 in Y
            SDR 003         // |di|*10^-3
            >REC
            cSTO .04        // cSTO ei

            // estimate derivative
            cRCL+ .02       // RCL+ zi. zi + ei
            TOP?            // show progress
              XEQ 91
            XEQ 84          // f(zi + ei)
              GTO 49        // if result is special or root, exit
            cRCL/ .06       // cRCL/ f(zi). f(zi + ei)/f(zi)
            cCHS            // -f(zi + ei)/f(zi)
            INC X           // 1 - f(zi + ei)/f(zi)
            cx=0?           // f(zi + ei) = f(zi)?
              GTO 49        // then exit (converged)
            cRCL .04        // cRCL ei
            cRCL/ Z         // ei/(1 - f(zi + ei)/f(zi))
            cRCL+ .02       // cRCL+ zi. zi + ei/(1 - f(zi + ei)/f(zi))
            cSTO .12        // cSTO wi. New wi

            // 2nd pass, Ostroswki's refine
            TOP?            // show progress
              XEQ 91
            XEQ 84          // f(wi)
              GTO 49        // if result is special or root, exit
            cRCL .06        // cRCL f(zi)
            cRCL/ Z         // f(zi)/f(wi)
            DEC X
            DEC X           // f(zi)/f(wi) - 2
            cRCL .12        // cRCL wi
            cRCL- .02       // cRCL- zi. wi - zi
            cRCL/ Z         // (wi - zi)/(f(zi)/f(wi) - 2)
            cRCL+ .12       // wi + (zi - wi)/(2 - f(zi)/f(wi)) (zi+1)
            SPEC?           // if f(zi) == 2*f(wi) omit 2nd pass by
              cRCL .12      // recalling wi

            // compute new |di| for next derivative estimation and check
            // convergence
            cFILL           // stack filled with zi+1
            cRCL- .02       // cRCL- zi. zi+1 - zi
            cABS            // |zi+1 - zi| (|d_{i+1}|)
            x=0?            // zi+1 = zi?
              GTO 49        // then exit (converged)
            cRCL .02        // cRCL zi
            cABS            // |zi|
            SDR 008         // |zi|*10^-8
            # 001
            SDR 008         // 10^-8
            MAX             // max(|zi|*10^-8, 10^-8)
            x>? Z           // |zi+1 - zi| < max(|zi|*10^-8, 10^-8)?
              SF .01        // yes, algorithm is converging
            x<> Z           // |zi+1 - zi|
            FC? .01         // converging?
              GTO 28        // no, skip convergence test
            x>? .01         // converged?
              GTO 49        // then exit (converged)
          LBL 28
            STO .01         // STO |di|

            // check user interrupt
            KEY? X          // any key?
              GTO 35        // no, continue
            # 035           // keycode for [<-]
            x=? Y           // key is [<-]?
              GTO 42        // yes, mark as no root found and exit
            cDROP

          LBL 35
            cRCL A          // cRCL zi+1
            // check iteration limit #
            DSZ .14         // DSZ counter
              GTO 14        // goto main loop
          LBL 42
            SF .00          // flag no root found

            // root found (F.00 reset) or not (F.00 set)
          LBL 49
            FS? .00         // no root?
              GTO 77        // no root, skip post process

            // beautify x.x...9999x...x and x.x...0000x...x roots
            # 012
            DBL?
              # 030
            SDR 003
            INC X
            STO .14         // ISG counter (1 to 12 (DBLOFF) or 30 (DBLON))
            # 003
            +
            STO .01         // di will always be counter + 3
            cRCL .02        // cRCL zi
            cSTO .10        // save zi in arg
          LBL 56            // beautifier loop
            RCL .10         // RCL arg
            RSD[->].14      // Re(zi) rounded to CNTR digits
            x<> .10
            RSD[->].01      // Re(zi) rounded to 3 more digits
            x!=? .10        // if both roundings are different then
              RCL L         // keep number intact (trying other rounding)
            STO .10         // update arg
            RCL .11         // RCL i_arg
            RSD[->].14      // Im(zi) rounded to CNTR digits
            x<> .11
            RSD[->].01      // Im(zi) rounded to 3 more digits
            x!=? .11        // if both roundings are different then
              RCL L         // keep number intact (but try other rounding)
            STO .11         // update i_arg
            INC .01         // increment pointers
            ISG .14
              GTO 56        // and loop

            // beautify almost real/imaginary roots
            FS? .03         // not needed in real mode
              GTO 70
            cRCL .10        // cRCL zi (partially beautified)
            EXPT
            x<> Y
            EXPT
            -               // orders of magnitude between Re and Im
            # 008
            x<? Y           // almost real?
              GTO 63        // yes, convert to real
            CHS
            x<? Y           // almost imaginary?
              GTO 70        // no, keep intact. Check new zi
            # 000           // here if almost imaginary, clear Re
            STO .10
            GTO 70          // zi beautified, check it
          LBL 63
            # 000           // here if almost real, clear Im
            STO .11

            // check if beautified zi is a better root
          LBL 70            // check beautified root
            cRCL .10        // cRCL zi (beautified)
            cx=? .02        // beautified zi = zi?
              GTO 77        // yes, no need to check, skip
            TOP?            // show progress
              XEQ 91
            XEQ 84          // compute f(zi) for beautified zi
              GTO 77        // beautified zi lets f(zi) special, skip...
                            // or is a true root, also skip
            cENTER
            cABS            // abs(f(zi)) with zi beautified
            cRCL .06
            cABS            // abs(f(zi)) with zi non-beautified
            x<? Z           // the beautified zi is no better?
              GTO 77        // yes, keep root intact
            cRCL A          // yes, root is the new zi. cRCL f(zi)
            cSTO .06        // cSTO f(zi)
            cRCL .10        // cRCL zi (beautified)
            cSTO .02        // cSTO zi (beautified)

            // build return stack
          LBL 77            // exit
            cRCL .08        // cRCL z0
            cSTO L          // input argument copied to LI (in real mode
                            // I is kept intact, as it was saved previously)
            CLSTK           // tidy stack (A-D lost)
            cRCL .06        // cRCL f(zi)
            cRCL .02        // cRCL zi. Now f(zi), zi on stack
            FC? .02         // restore user's D flag
              CF D
            RCLM .00        // RCLM mode, restore user mode
            FS? .03         // adjust stack for the real case?
              <> XZYY       // yes, adjust...
            TOP?
              XEQ 91        // show result
            FC? .00         // root found?
              RTN           // yes, that's enough
            TOP?
              MSG 20        // no, report it
          RTN+1             // to behave like a test

            // wrapper for program to be solved. Saves the argument,
            // fills the stack and checks return value. Behaves like
            // a test. The test returns false if the return value is
            // neither special nor zero. Returns true otherwise,
            // setting flag .00 if return value is special or storing
            // 0 + 0j in f(zi) and the argument in zi if the argument is
            // a root
          LBL 84
            cSTO .10        // cSTO arg. Save argument
            cFILL
            FS? .03         // real mode?
              FILL
            FC? .02         // honour user's D flag
              CF D
            XEQa            // evaluate f
            SF D            // restore D flag
            FS? .03         // force 0 imaginary part in real mode
              # 000
            FS? .03
              x<> Y
            SF .00
            SPEC?           // Re(f) is special?
              RTN           // yes, test returns true & F.00 set
            x<> Y
            SPEC?           // Im(f) is special?
              RTN           // yes, test returns true & F.00 set
            x<> Y
            CF .00
            cx!=0?          // argument is not a root?
              RTN+1         // not a root, test returns false
            cSTO .06        // cSTO f(zi). Is a root, save 0+0j in f(zi)
            cRCL .10        // cRCL arg
            cSTO .02        // cSTO zi. Store the root in zi
            RTN             // test returns true & F.00 reset

          LBL 91            // show progress
            FS? .03
              VIEW X
            FC? .03
              cVIEW X
            RTN

        END

César - Information must flow.
Find all posts by this user
Quote this message in a reply
04-16-2017, 05:30 PM
Post: #2
RE: (WP-34S) Real/complex solver
And this is the v1.6L:

Code:

/*
* RCSV, an Ostrowski-Traub based real/complex solver
*
*   v1.6L (20170414) by M. César Rodríguez GGPL
*
* Usage:
*   When used as a real solver (called as 'RSV'):
*     Input:
*       X:      contains an initial root estimation
*       alpha:  contains the name of the program to be solved. This program
*               is called with the stack (an SSIZE8 one) filled with an
*               estimation of the root. The program must return a number in X
*     Output:
*       X:      the found root or the value of the last estimation for the root
*       Y:      the value of the program to be solved evaluated at the returned
*               value in X
*       ZT:     lost (zeroed)
*       L:      the input X value
*       IJK:    kept intact
*       ABCD:   lost (zeroed)
*
*   When used as a complex solver (called as 'CSV'):
*     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
*       JK:     kept intact
*       ABCD:   lost (set to zero)
*
*     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 (400 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. The same will occur if the
*     program to be solved throws an error (e.g. trying to access a non
*     existent register)
*
*   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 for all internal calculations in order not to throw
*       errors when a calculation returns Infinite or NaN, but the program
*       to be solved is called with the flag D the user set. The D flag
*       is also restored at 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 approximated 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^-8, 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 (0 in real mode)
*     - 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 400 iterations
*     - the user interrupts the program pressing [<-]
*
*   When a root is found, it is "beautified" by:
*     - changing numbers like 1.124999973587 or 1.125000032518
*       to 1.125. The code looks for chains of four or more 0s or 9s and
*       rounds the number accordingly
*     - in complex mode, 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, in complex mode, 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
*
* 217 steps, 15 local registers
*
*/
          LBL'RSV'          // Real SolVer
            LocR 15
            SF .03          // flag the real mode
            JMP CSV_entry

          LBL'CSV'          // Complex SolVer
            LocR 15         // Local registers and flags:
                            //  .00       mode bits
                            //  .01       |d_i| (that's |z_i - z_{i-1}|)
                            //  .02-03    z_i
                            //  .04-05    e_i
                            //  .06-07    f(z_i)
                            //  .08-09    z_0
                            //  .10-11    arg, argument for program to
                            //            solve, also beautified zi
                            //  .12-13    w_i
                            //  .14       iteration counter
                            //  F.00      no convergence, or f evaluates
                            //            to special, or user interrupt
                            //            ([<-])
                            //  F.01      algorithm is converging
                            //  F.02      backup of D flag
                            //  F.03      true in real mode

CSV_entry:: STOM .00        // STOM mode. Mode is restored at exit
            SSIZE8
            FS? D           // backup flag D
              SF .02
            SF D            // specials are not errors here
            FS? .03         // in real mode, I must be maintained, so
              y<> I         // it is saved as the imaginary part of Z_0
            cSTO .08        // cSTO z0. Initial guess
            FS? .03         // in real mode, the imaginary part of the
              # 000         // estimates should be 0
            FS? .03
              x<> Y
            cSTO .02        // cSTO zi. Estimate
            cABS            // |z0|
            x=0?            // |z0| = 0?
              # 001
            SDR 003         // |d0| = |z0|*10-3, or 10^-3 if |z0| = 0
            STO .01         // later it is divided again by 10^3. STO |di|
            # 040
            SDL 001         // 400, maximum iter #
            STO .14         // STO counter
            cRCL .02        // cRCL zi

            // first pass, Newton's estimate (zi on stack)
CSV_mloop:: cSTO .02        // cSTO zi
            TOP?            // show progress
              GSB CSV_msg
            GSB CSV_eval    // eval function, check for special or root
              JMP CSV_exit  // if result is special or root, exit
            cSTO .06        // cSTO f(zi)

            // compute e_i to estimate derivative
            RCL .03         // RCL Im(zi)
            # 000
            x=? Y           // zi is real?
              JMP CSV_noai  // yes, e_i must also be real
            ACOS            // no. pi/2 (or 90º or 100 grad)
            # 004
            *
            RAN#
            *               // ai, random angle between 0 and 2pi
CSV_noai::  RCL .01         // RCL |di|. ai or 0 in Y
            SDR 003         // |di|*10^-3
            >REC
            cSTO .04        // cSTO ei

            // estimate derivative
            cRCL+ .02       // RCL+ zi. zi + ei
            TOP?            // show progress
              GSB CSV_msg
            GSB CSV_eval    // f(zi + ei)
              JMP CSV_exit  // if result is special or root, exit
            cRCL/ .06       // cRCL/ f(zi). f(zi + ei)/f(zi)
            cCHS            // -f(zi + ei)/f(zi)
            INC X           // 1 - f(zi + ei)/f(zi)
            cx=0?           // f(zi + ei) = f(zi)?
              JMP CSV_exit  // then exit (converged)
            cRCL .04        // cRCL ei
            cRCL/ Z         // ei/(1 - f(zi + ei)/f(zi))
            cRCL+ .02       // cRCL+ zi. zi + ei/(1 - f(zi + ei)/f(zi))
            cSTO .12        // cSTO wi. New wi

            // 2nd pass, Ostroswki's refine
            TOP?            // show progress
              GSB CSV_msg
            GSB CSV_eval    // f(wi)
              JMP CSV_exit  // if result is special or root, exit
            cRCL .06        // cRCL f(zi)
            cRCL/ Z         // f(zi)/f(wi)
            DEC X
            DEC X           // f(zi)/f(wi) - 2
            cRCL .12        // cRCL wi
            cRCL- .02       // cRCL- zi. wi - zi
            cRCL/ Z         // (wi - zi)/(f(zi)/f(wi) - 2)
            cRCL+ .12       // wi + (zi - wi)/(2 - f(zi)/f(wi)) (zi+1)
            SPEC?           // if f(zi) == 2*f(wi) omit 2nd pass by
              cRCL .12      // recalling wi

            // compute new |di| for next derivative estimation and check
            // convergence
            cFILL           // stack filled with zi+1
            cRCL- .02       // cRCL- zi. zi+1 - zi
            cABS            // |zi+1 - zi| (|d_{i+1}|)
            x=0?            // zi+1 = zi?
              JMP CSV_exit  // then exit (converged)
            cRCL .02        // cRCL zi
            cABS            // |zi|
            SDR 008         // |zi|*10^-8
            # 001
            SDR 008         // 10^-8
            MAX             // max(|zi|*10^-8, 10^-8)
            x>? Z           // |zi+1 - zi| < max(|zi|*10^-8, 10^-8)?
              SF .01        // yes, algorithm is converging
            x<> Z           // |zi+1 - zi|
            FC? .01         // converging?
              JMP CSV_updt  // no, skip convergence test
            x>? .01         // converged?
              JMP CSV_exit  // then exit (converged)
CSV_updt::  STO .01         // STO |di|

            // check user interrupt
            KEY? X          // any key?
              JMP CSV_nok   // no, continue
            # 035           // keycode for [<-]
            x=? Y           // key is [<-]?
              JMP CSV_flag  // yes, mark as no root found and exit
            cDROP

CSV_nok::   cRCL A          // cRCL zi+1
            // check iteration limit #
            DSZ .14         // DSZ counter
            JMP CSV_mloop   // goto main loop
CSV_flag::  SF .00          // flag no root found

            // root found (F.00 reset) or not (F.00 set)
CSV_exit::  FS? .00         // no root?
              JMP CSV_nopp  // no root, skip post process

            // beautify x.x...9999x...x and x.x...0000x...x roots
            # 012
            DBL?
              # 030
            SDR 003
            INC X
            STO .14         // ISG counter (1 to 12 (DBLOFF) or 30 (DBLON))
            # 003
            +
            STO .01         // di will always be counter + 3
            cRCL .02        // cRCL zi
            cSTO .10        // save zi in arg
CSV_bloop:: RCL .10         // RCL arg
            RSD[->].14      // Re(zi) rounded to CNTR digits
            x<> .10
            RSD[->].01      // Re(zi) rounded to 3 more digits
            x!=? .10        // if both roundings are different then
              RCL L         // keep number intact (trying other rounding)
            STO .10         // update arg
            RCL .11         // RCL i_arg
            RSD[->].14      // Im(zi) rounded to CNTR digits
            x<> .11
            RSD[->].01      // Im(zi) rounded to 3 more digits
            x!=? .11        // if both roundings are different then
              RCL L         // keep number intact (but try other rounding)
            STO .11         // update i_arg
            INC .01         // increment pointers
            ISG .14
              JMP CSV_bloop // and loop

            // beautify almost real/imaginary roots
            FS? .03         // not needed in real mode
              JMP CSV_chkb
            cRCL .10        // cRCL zi (partially beautified)
            EXPT
            x<> Y
            EXPT
            -               // orders of magnitude between Re and Im
            # 008
            x<? Y           // almost real?
              JMP CSV_real  // yes, convert to real
            CHS
            x<? Y           // almost imaginary?
              JMP CSV_chkb  // no, keep intact. Check new zi
            # 000           // here if almost imaginary, clear Re
            STO .10
            JMP CSV_chkb    // zi beautified, check it
CSV_real::  # 000           // here if almost real, clear Im
            STO .11

            // check if beautified zi is a better root
CSV_chkb::  cRCL .10        // cRCL zi (beautified)
            cx=? .02        // beautified zi = zi?
              JMP CSV_nopp  // yes, no need to check, skip
            TOP?            // show progress
              GSB CSV_msg
            GSB CSV_eval    // compute f(zi) for beautified zi
              JMP CSV_nopp  // beautified zi lets f(zi) special, skip...
                            // or is a true root, also skip
            cENTER
            cABS            // abs(f(zi)) with zi beautified
            cRCL .06
            cABS            // abs(f(zi)) with zi non-beautified
            x<? Z           // the beautified zi is no better?
              JMP CSV_nopp  // yes, keep root intact
            cRCL A          // yes, root is the new zi. cRCL f(zi)
            cSTO .06        // cSTO f(zi)
            cRCL .10        // cRCL zi (beautified)
            cSTO .02        // cSTO zi (beautified)

            // build return stack
CSV_nopp::  cRCL .08        // cRCL z0
            cSTO L          // input argument copied to LI (in real mode
                            // I is kept intact, as it was saved previously)
            CLSTK           // tidy stack (A-D lost)
            cRCL .06        // cRCL f(zi)
            cRCL .02        // cRCL zi. Now f(zi), zi on stack
            FC? .02         // restore user's D flag
              CF D
            RCLM .00        // RCLM mode, restore user mode
            FS? .03         // adjust stack for the real case?
              <> XZYY       // yes, adjust...
            TOP?
              GSB CSV_msg   // show result
            FC? .00         // root found?
              RTN           // yes, that's enough
            TOP?
              MSG 20          // no, report it
          RTN+1             // to behave like a test

            // wrapper for program to be solved. Saves the argument,
            // fills the stack and checks return value. Behaves like
            // a test. The test returns false if the return value is
            // neither special nor zero. Returns true otherwise,
            // setting flag .00 if return value is special or storing
            // 0 + 0j in f(zi) and the argument in zi if the argument is
            // a root
CSV_eval::  cSTO .10        // cSTO arg. Save argument
            cFILL
            FS? .03         // real mode?
              FILL
            FC? .02         // honour user's D flag
              CF D
            XEQa            // evaluate f
            SF D            // restore D flag
            FS? .03         // force 0 imaginary part in real mode
              # 000
            FS? .03
              x<> Y
            SF .00
            SPEC?           // Re(f) is special?
              RTN           // yes, test returns true & F.00 set
            x<> Y
            SPEC?           // Im(f) is special?
              RTN           // yes, test returns true & F.00 set
            x<> Y
            CF .00
            cx!=0?          // argument is not a root?
              RTN+1         // not a root, test returns false
            cSTO .06        // cSTO f(zi). Is a root, save 0+0j in f(zi)
            cRCL .10        // cRCL arg
            cSTO .02        // cSTO zi. Store the root in zi
            RTN             // test returns true & F.00 reset

CSV_msg::   FS? .03         // show progress
              VIEW X
            FC? .03
              cVIEW X
            RTN

        END

César - Information must flow.
Find all posts by this user
Quote this message in a reply
Post Reply 




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