(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
|
Search
Member List
Calendar
Help
