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