[34S] Time diff. STO.00 vs. STO 00 ?

[Sikuq]

(03-13-2014 09:53 PM)Paul Dale Wrote:  

(03-13-2014 09:20 PM)Sikuq Wrote:  However I always wondered just how the 34S actually calculates say a square root or the Gamma function. Does it use a power series, other basic functions, or say Fourier? And did you get those functions from somewhere else and just compiled them in there?

Square root, the four arithmetic operations, natural logarithm and exponential are included in the decNumber library we used. I replaced the natural logarithm (twice) because the supplied version was unacceptably slow, I'm still not entirely happy with its performance. For reference, square root uses Newton's method after an initial guess -- this is a fairly standard technique.

Pretty much everything else I wrote after finding suitable standard algorithms or in a couple of cases doing my own expansions. As an example of the former, the gamma function uses a rapidly converging series -- unfortunately I don't remember the source I used, but I will have the original paper stored away somewhere. I ended up purchasing a not-insignificant number of books about numeric analysis and transcendental function implementation and spent lots of time searching the Internet for suitable algorithms.

If you look in doc/formulas there are some erratic notes about what I used where and doc/distribution-formulas contains a longish exchange between Dieter and others giving incremental improvements to the statistical distributions. These documents were primarily for my use and aren't well organised but I figured some record was better than none.

Finally, a few of the functions have been examined in detail by members of this forum and sometimes, changes and improvements have been added. For example, Dieter's work on getting the Lambert W function accurate.


- Pauli

Pauli,

Despite Walter's possible consternation I allow myself to re-quote you fully. Your post is so profoundly interesting to me and I am sure most other 34S users that it would be a shame to truncate it. I wish I had had the knowledge and involvement at that time to experience what I can only describe and imagine as a transcontinental total brainstorming which at times lead to Nobel Prizes. I can vividly imagine your discussions as to size, accuracy, speed, convergence, etc. I never did know that you guys were so involved at that fundamental level. My interest in this topic has boosted and I must look into your sources and elsewhere.

Years ago I did my own simple Newton-Raphson and secant method solve routine on my HP41C and I remember it was sometimes calculating overnight before it printed out on my small HP thermal printer. I did have two of the HP 41C IL-Loop tape drives but I can't remember if I involved those. I was solving some huge complicated Gamma functions. One of my problems was to discern which of the multiple solutions were legitimate. I still have problems with those issues; I can never make up my mind as I have often been unable to establish a delimiting criteria for solution detection. Thus I often worry about the consequences of the solution I let go of coming home to roost.

I'm going to have a sleepless night thinking about your exciting endeavors.

- Sikuq

[walter b]

Dear Sikuq,

You're right guessing is another hobby of me. But now - since it looks like you're a lonely offspring of a (not as lonely!) medieval Viking monk on a lonely island I can only congratulate to your location: Enjoy and keep it! ((Just some reasoning: regards copying as a value per se increasing the value of the copied text (compare e.g. Umberto Eco's "The Name of the Rose"), long distance given in nmi to a university, chosen forum name (sounds a bit Inuktitut to me), lack of grammatical knowledge (contact to Latin lost for long), etc. ))

Seriously, I didn't want to offend you - your questions just sound sometimes that a little bit of own research would ... ummh ... help. But times have changed and are changing - nowadays the intelligence is in the web ...

d;-)

I'll try to resist to comment on your future posts but can't promise to succeed

[Paul Dale]

The majority of the discussions about the numerical accuracy are stored in the old forum archives. There was surprisingly little discussion about this elsewhere.


- Pauli

[Sikuq]

Walter,

I am not one bit offended, au contraire! I am in fact honored by the discussions going on here. Please keep commenting. It does force me to do some new research and the feedback by Marcus, yourself, Pauli, Franz, and others are always very informative. Your WP34S team is such an unusual group of such high intellect that at times I feel like a fool.

You are this time amazingly correct in you assumptions except for the lonely offspring part. I did go west like Vikings do and I did snatch a cute little local girl like Vikings do...

- Sikuq

Edit 1,2.

Walter, f.y.i., Sikuq= ice in my Eskimo variant. You are very observant indeed.

[Sikuq]

(03-14-2014 07:18 AM)Paul Dale Wrote:  The majority of the discussions about the numerical accuracy are stored in the old forum archives. There was surprisingly little discussion about this elsewhere.


- Pauli

Pauli,

The story behind your 34S is astounding; yes, it certainly has not been exposed to the extent I would have guessed.

I woke up several times over night thinking of the 34S team running around in Australia and Germany seeking algorithm books, and also recalling my own attempts to make my own Newton routine faster. Sometimes I would simple-graph the function on the HP41C loop printer and count the roots.

Does the user FIX setting in any way affect the actual step length used?

- Sikuq

PS: Pauli, is there any way you can test how many loops it takes to calculate the square-root of say 1234567 on the 34S? Is the step length fixed or does it reduce itself closer to convergence (i.e. sign change)? And how did it calculate the initial guess?

[Thomas Klemm]

(03-14-2014 05:01 PM)Sikuq Wrote:  PS: Pauli, is there any way you can test how many loops it takes to calculate the square-root of say 1234567 on the 34S? Is the step length fixed or does it reduce itself closer to convergence (i.e. sign change)? And how did it calculate the initial guess?

You might want to have a look at how to calculate the fast inverse square root. But I must admit that I have no idea whether this is used in the WP-34s.
The algorithm used in most HP calculators is similar to the digit-by-digit method you might have learned at school but uses only addition and shifts. Thus just counting the amount of loops isn't enough if you are interested in performance.

Kind regards
Thomas

[Thomas Klemm]

(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                                                           */
/* ------------------------------------------------------------------ */

[Sikuq]

(03-14-2014 05:34 PM)Thomas Klemm Wrote:  

(03-14-2014 05:01 PM)Sikuq Wrote:  PS: Pauli, is there any way you can test how many loops it takes to calculate the square-root of say 1234567 on the 34S? Is the step length fixed or does it reduce itself closer to convergence (i.e. sign change)? And how did it calculate the initial guess?

You might want to have a look at how to calculate the fast inverse square root. But I must admit that I have no idea whether this is used in the WP-34s.
The algorithm used in most HP calculators is similar to the digit-by-digit method you might have learned at school but uses only addition and shifts. Thus just counting the amount of loops isn't enough if you are interested in performance.

Kind regards
Thomas

Thomas,

I thought I had posted a response to you earlier but it appears to have been previewed only - sorry.

Thank you for your two posts and the code which I am studying. The inverse function was entirely new to me.

I am curious as to just why the loop-cycles are not a measure of efficiency given the same accuracy? As I wrote to Pauli I am astounded to learn that the 34S and other computers use the standard Newton tangent way to solve roots, a way I understand is prominent.

Zurich I read was voted the most livable city, lucky you.

Thanks.

- Sikuq

[Sikuq]

(03-14-2014 07:18 AM)Paul Dale Wrote:  The majority of the discussions about the numerical accuracy are stored in the old forum archives. There was surprisingly little discussion about this elsewhere.
- Pauli

Pauli,

I've looked through some of your earlier posts in regards to algorithms. That led me to the square root algorithm by William E. Egbert from HP while working on his HP67 project as described in the HP Journal issue from May 1977. Is that the one you use in the WP34S or is there a better one? It certainly looks like a fast way of getting to the result using only basic WP34S functions.

- Sikuq

[Thomas Klemm]

(03-17-2014 05:12 PM)Sikuq Wrote:  Is that the one you use in the WP34S or is there a better one?

They are not the same: Newton vs. digit-by-digit

Code:

/* decNumberSquareRoot -- square root operator                        */

/* 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                                                           */

Cheers
Thomas

[Sikuq]

(03-17-2014 06:14 PM)Thomas Klemm Wrote:  

(03-17-2014 05:12 PM)Sikuq Wrote:  Is that the one you use in the WP34S or is there a better one?

They are not the same: Newton vs. digit-by-digit

Code:

/* decNumberSquareRoot -- square root operator                        */

/* 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                                                           */

Cheers
Thomas

Thomas,

Many thanks. This algorithm stuff gets huge quickly as size is at a premium on WP34S I'd suppose. I'm tempted to key-in the Egbert routines and check convergence times versus the Newton as indicated by you although I have a feeling I would be at a disadvantage using program lines versus pre-compiled c-code? Am I correct?

- Sikuq

[Thomas Klemm]

(03-14-2014 09:09 PM)Sikuq Wrote:  I am curious as to just why the loop-cycles are not a measure of efficiency given the same accuracy?

In case of Newton's method you may not need many iterations but time consuming operations like a division are used. The digit-by-digit method may use more cycles but the operations are simple and
fast. Thus it really depends on what is going on within that loop.

HTH
Thomas

[Sikuq]

(03-17-2014 07:38 PM)Thomas Klemm Wrote:  

(03-14-2014 09:09 PM)Sikuq Wrote:  I am curious as to just why the loop-cycles are not a measure of efficiency given the same accuracy?

In case of Newton's method you may not need many iterations but time consuming operations like a division are used. The digit-by-digit method may use more cycles but the operations are simple and
fast. Thus it really depends on what is going on within that loop.

HTH
Thomas

Thomas,

Thanks for the heads-up. It makes me appreciate just how much work they must have done when they stuffed all those algorithms into the WP34S.