(03-14-2014 05:34 PM)Thomas Klemm Wrote: But I must admit that I have no idea whether this is used in the WP-34s.
Meanwhile I had a look at the source code. This comment is from the file
wp34s-code/decNumber/decNumber.c:
Code:
/* ------------------------------------------------------------------ */
/* decNumberSquareRoot -- square root operator */
/* */
/* This computes C = squareroot(A) */
/* */
/* res is C, the result. C may be A */
/* rhs is A */
/* set is the context; note that rounding mode has no effect */
/* */
/* C must have space for set->digits digits. */
/* ------------------------------------------------------------------ */
/* This uses the following varying-precision algorithm in: */
/* */
/* Properly Rounded Variable Precision Square Root, T. E. Hull and */
/* A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
/* pp229-237, ACM, September 1985. */
/* */
/* The square-root is calculated using Newton's method, after which */
/* a check is made to ensure the result is correctly rounded. */
/* */
/* % [Reformatted original Numerical Turing source code follows.] */
/* function sqrt(x : real) : real */
/* % sqrt(x) returns the properly rounded approximation to the square */
/* % root of x, in the precision of the calling environment, or it */
/* % fails if x < 0. */
/* % t e hull and a abrham, august, 1984 */
/* if x <= 0 then */
/* if x < 0 then */
/* assert false */
/* else */
/* result 0 */
/* end if */
/* end if */
/* var f := setexp(x, 0) % fraction part of x [0.1 <= x < 1] */
/* var e := getexp(x) % exponent part of x */
/* var approx : real */
/* if e mod 2 = 0 then */
/* approx := .259 + .819 * f % approx to root of f */
/* else */
/* f := f/l0 % adjustments */
/* e := e + 1 % for odd */
/* approx := .0819 + 2.59 * f % exponent */
/* end if */
/* */
/* var p:= 3 */
/* const maxp := currentprecision + 2 */
/* loop */
/* p := min(2*p - 2, maxp) % p = 4,6,10, . . . , maxp */
/* precision p */
/* approx := .5 * (approx + f/approx) */
/* exit when p = maxp */
/* end loop */
/* */
/* % approx is now within 1 ulp of the properly rounded square root */
/* % of f; to ensure proper rounding, compare squares of (approx - */
/* % l/2 ulp) and (approx + l/2 ulp) with f. */
/* p := currentprecision */
/* begin */
/* precision p + 2 */
/* const approxsubhalf := approx - setexp(.5, -p) */
/* if mulru(approxsubhalf, approxsubhalf) > f then */
/* approx := approx - setexp(.l, -p + 1) */
/* else */
/* const approxaddhalf := approx + setexp(.5, -p) */
/* if mulrd(approxaddhalf, approxaddhalf) < f then */
/* approx := approx + setexp(.l, -p + 1) */
/* end if */
/* end if */
/* end */
/* result setexp(approx, e div 2) % fix exponent */
/* end sqrt */
/* ------------------------------------------------------------------ */
Search
Member List
Calendar
Help
/
[34S] Time diff. STO.00 vs. STO 00 ?
