The Museum of HP Calculators

HP Forum Archive 06

[ Return to Index | Top of Index ]

More trig--specifically speed vs accuracy
Message #1 Posted by Cameron on 22 Oct 2001, 11:51 a.m.

All the recent discussions about trig functions on calculators that don't support them natively have got my inquiry juices flowing. I devoured Viktor Toth's recent contribution for the 12C with interest. The timings in particular caught my eye. (Thanks for posting that stuff in the articles forum Viktor).

I don't have a clue what algorithm Viktor is using and I'm not familiar with the 12C. However that didn't deter me from implementing a Taylor Polynomial algorithm on my 16C. The results were both disappointing and revealing.

I managed sine in 94 instructions. To put that in perspective, the sine routine uses 37 instructions. The rest are taken up with computing y**x and x!. I could probably shave two handfulls of instructions off if I took out the error detection and optimised register usage. Even so, it would still be a weighty program.

The most disappointing thing is that it takes around 155 seconds to produce a result. Which leads me to the point of this posting. What are typical times for sine on those Voyagers (10C, 11C and 15C) that support the function natively? Since my 1980's Casio algebraic does it in less than 2 seconds, I assume the Voyagers are in the same ball park. Assuming that they are, what algorithm do those Voyagers use?

I'm aware of CORDIC but I can't see a way of implementing it (easily) on the 16C so I can do the comparison, largely because the bit shift operators are not available in floating point mode.

I'll try to reduce the time of my implementation by precomputing the Taylor coefficients. The draw-back to that approach is that it requires 8 * 7 bytes of register storage. I'll post comparative timings when I've done that.

Your thoughts?

Cameron

PS: my 16C simulator does the 17th degree Taylor Polynomial algorithm in less than the key debounce time on my stopwatch. That's on a 750 MHz processor. What was the processor clock speed of the 16C? 8 KHz?

      
Re: More trig--specifically speed vs accuracy
Message #2 Posted by W. Bruce Maguire II on 22 Oct 2001, 12:15 p.m.,
in response to message #1 by Cameron

To answer some questions,
Victor: could you post something to let us know what algorithm(s) you are using?

As to the timing question, I can only add a little information. I'm not sure that "modern" calculators still use the CORDIC algorithm---they may be using something a little faster. But, when the functions are native to the calculator, they are in firmware---not "interpreted" line-by-line as in the key-programmable calculators. My (only slightly-informed) guess is that this could easily account for TWO orders of magnitude in speed.

My 1.5 cents,
Bruce.

            
Re: More trig--specifically speed vs accuracy
Message #3 Posted by Viktor Toth on 22 Oct 2001, 9:39 p.m.,
in response to message #2 by W. Bruce Maguire II

Bruce,

I used the Taylor expansions, employing some of the tricks Ex-PPC member mentioned in order to avoid unnecessary calculations.

Viktor

                  
Not to mention those handy STO instructions
Message #4 Posted by Cameron on 23 Oct 2001, 1:47 a.m.,
in response to message #3 by Viktor Toth

Those STO <op> <regnum> instructions are a real bonus!

Am I correct in assuming that they apply the operator to X and the current content of regnum and then store the result back in regnum? Something like:

RCL <regnum>
<op>
STO <regnum>

If so, they are so elegant; so succinct. I can't see how they'd take up much ROM space. I wonder why the 16C doesn't support them? There are 11 instances of them in Viktor's program, which would expand to 33 instructions on the 16C.

Cameron

                        
Re: Not to mention those handy STO instructions
Message #5 Posted by Viktor Toth on 23 Oct 2001, 7:30 a.m.,
in response to message #4 by Cameron

Cameron,

Yes, register arithmetic IS a real bonus.

And your understanding of their function is largely correct, but there is an important point you might have missed. STO+ 00 adds the contents of X to register 00 _without_ losing the contents of either T or LASTx.

Viktor

                              
Questions about register arithmetic.
Message #6 Posted by Cameron on 24 Oct 2001, 2:02 p.m.,
in response to message #5 by Viktor Toth

Thanks for the explanation, Viktor.

STO+ 00 adds the contents of X to register 00 _without_ losing the contents of either T or LASTx.

Is that a property of all operator/register combinations or only a special property of using register zero?

Also, does the operator extend beyond the four basic arithmetic ones? For all registers?

Cameron

                                    
Re: Questions about register arithmetic.
Message #7 Posted by Vieira, Luiz C. (Brazil) on 24 Oct 2001, 5:04 p.m.,
in response to message #6 by Cameron

Hello, Cameron.

Register arithmetic is available in almost all HewPack calcs. You asked about the operations: the four basic ones are the first target. By looking at the HP Museum Data Sheets for some calculators, a STO ^ is mentioned, but I did not dare finding any calculator that would raise to the X-register power the contents of a register.

About all registers access: if we consider the number of bytes per instruction (if they can be part of a program) then all register number should be considered. In the HP41C, arithmetic is available for the stack registers (ST<op> <X Y Z T L>), direct <00> to <99> and indirect to all mentioned (say, ST<op> IND <X Y Z T L 00 TO 99>. o.k.?

It grows better when we get to the HP15C. Arithmetic on registers accepts STO <op> <0 to 9 .0 to .9 I (i)> and RCL <op> <0 to 9 .0 to .9 I (i)>. That is, the contents of the X-register is the result of the previous contents <op> the register contents. Previous X-contents? to LASTx. Good, ahn?

The RCL <op> is available in the HP42S, too.

There is a bit more to say, but I have no manuals in hand, so before I write some incorrect info, I will wait for your other questions, right?

Hope it helps. Regards.

      
Re: More trig--specifically speed vs accuracy
Message #8 Posted by Ex-PPC member on 22 Oct 2001, 2:18 p.m.,
in response to message #1 by Cameron

Hi, Cameron:

You wrote:

"I managed sine in 94 instructions. To put that in perspective, the sine routine uses 37 instructions. The rest are taken up with computing y**x and x!. I could probably shave two handfulls of instructions off if I took out the error detection and optimised register Even so, it would still be a weighty program."

You don't need to emulate y^x nor x! or your HP-16C. The Taylor Series Expansion for sin(x) is:

sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...

If you compute this as a sum of N terms, the terms T(n) are:

T(1) = +x, T(2) = -x^3/3!, T(3) = +x^5/5!, T(4) = -x^7/7!, etc.

and it seems necessary to have both an y^x and an x! routines. But you can compute each term T(n) as a function of the previous one, using this simple relation:

T(n+1) = T(n) * [ -x*x/((2n-1)*(2n-2)) ]

So, for instance: T(4) = T(3)*[ -x*x/(7*6)] = -x^5/5 * (-x^2/(6*7)) = +x^7/7!

as it should. Using this relation, you can simply keep the previous term T(n) somewhere, be it the stack or a register, then compute the next term, T(n+1) by simply multiplying it times -x*x [which you'd do well to compute at first and store somewhere to avoid recomputing it in each iteration], then dividing by 2n-1 and then again by 2n-2. Once you've got this new term, T(n+1), add it to the ongoing sum, and keep it as the new term from which to compute the next.

That way you don't need y^x nor x!, and your program will much faster, and even shorter. With a little persevering, you should get a 20 step sine routine easily, and even shorter is possible. :-)

Hope this helps ...

            
Re: More trig--specifically speed vs accuracy
Message #9 Posted by Andrés C. Rodríguez (Argentina) on 23 Oct 2001, 9:28 a.m.,
in response to message #8 by Ex-PPC member

Mr. Ex-PPC member:

Thank you for an interesting posting; good reading even for people like me who are not particularly interested in the "trigs" challenge at this time.

                  
Re: More trig--specifically speed vs accuracy
Message #10 Posted by Ex-PPC member on 23 Oct 2001, 10:02 a.m.,
in response to message #9 by Andrés C. Rodríguez (Argentina)

Mr. Andres C. Rodriguez wrote: "Thank you for an interesting posting; good reading even for people like me who are not particularly interested in the "trigs" challenge at this time."

You're welcome. I'm almost finished with my version, just optimizing it a little and trying to include extra functionality in the free space left, but the issue is essentially solved as far as I'm concerned and I'm moving on to other interesting hp-calc programming topics. Thanks for your always balanced and kind comments.

            
Thanks and <blush> help...
Message #11 Posted by Cameron on 24 Oct 2001, 11:13 a.m.,
in response to message #8 by Ex-PPC member

Thanks Ex-PPC. I've not yet risen to the challenge due to irritating real-world intrusions. Could you put an old bloke, who's brain has gone to mush (not from CJD but from working with Microsoft products) out of his misery? What's the equivalent relation for cosine. I'm embarrassed to confess that the first principles work is beyond my feeble intellect.

TIA

Cameron

                  
Re: Thanks and <blush> help...
Message #12 Posted by Ex-PPC member on 24 Oct 2001, 11:35 a.m.,
in response to message #11 by Cameron

Hi again, Cameron. The Taylor series expansion for the cosine is:

cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ...

and the relationship between the term T(n) and the term T(n+1) is exactly the same as that for the sine terms ! Matter of fact the only difference is that T(1) is "1" for the cosine and "x" for the sine.

However, I would suggest that, once you've computed the sine, you simply compute the cosine using this well-known equation:

cos(x) = square_root(1-(sin(x))^2)

Hope that helps.

                        
Duh! <more red faces>
Message #13 Posted by Cameron on 24 Oct 2001, 11:47 a.m.,
in response to message #12 by Ex-PPC member

I knew the Taylor series for cosine. It was the T(n) to T(n+1) relationship I imagined I'd have a hard time deriving. I'm not going compound my public act of doltery and ask why. I will ponder it some.

As for the sine-cosine relation, ummm, I knew that there was one but couldn't dredge it up.

Thanks for taking the time to spell it out.

Cameron

                              
Re: <more red faces> - Not only you
Message #14 Posted by Vieira, Luiz C. (Brazil) on 24 Oct 2001, 12:31 p.m.,
in response to message #13 by Cameron

Hey, Cameron;

I believe you have been courageous enough expressing your doubts about all this strong demonstration of practice, usage of mathematics tools. I had um own doubts, but I fake dead...

I felt brought back to high-school, and felt myself good at least understanding it all. Building it all by myself? I was delving into a way to easily obtain the minus sign (-) for the odd terms, and the solution is far from these concerns... Amazing.

That’s what I call a multi-purpose forum. I’ve been in silence till now cause both a lack in spare time (just reading) and knowledge (just understanding has been enough for now...)

And... (may I borrow this space?) Mr. Ex-PPC Member, what other challenges are you chasing for now? You keep teasing us by just mentioning you are working on new programs... Please ??? Will you share with us?

Cheers.

                                    
Re: <more red faces> - Not only you
Message #15 Posted by Ex-PPC member on 24 Oct 2001, 1:21 p.m.,
in response to message #14 by Vieira, Luiz C. (Brazil)

Mr. Luiz wrote:

"And... (may I borrow this space?) Mr. Ex-PPC Member, what other challenges are you chasing for now? You keep teasing us by just mentioning you are working on new programs... Please ??? Will you share with us?"

Certainly. As my 99-step version of the trigonometric functions for the 12C will also run on the 16C almost unchanged (as I do not use storage arithmetic at all, STO+, STO-, etc), I'm now considering implementing on the HP-16C two important floating-point functions it does lack, namely y^x and square root, for arbitrary floating point arguments of course. Unlike the HP-12C, the HP-16C doesn't have either e^x or ln(x), so implementing y^x accurately and fast is quite a challenge indeed.

                                          
Allow me to return the favour...
Message #16 Posted by Cameron on 24 Oct 2001, 1:45 p.m.,
in response to message #15 by Ex-PPC member

...and share my modest attempt at y**x on the 16C.

A Subroutine to Compute Y**X

I considered doing X**Y but that doesn't seem to be the HP "way".

LBL 2 CF 0 ; Prep flag 0 for use. x == 0 ; X == 0? SF 0 ; then set flag zero [1] F? 0 RTN

RDN ; move the 1 out of the way

x > 0 ; +ve power? (we've already handled x == 0) GTO C

RUP ; get that 1 2 DIV x == 0 RTN ; cos we're in integer mode

CLx ; disable lift [1] RDN CHS RDN 1/X RUP

LBL C

RUP ; return 1 to X SWAP XY ; power into X; 1 into Y STO I ; power is iterator RDN ; discard power SWAP XY ; base into X; 1 into Y MUL ; discard 1; establish X as constant

LBL D ; start loop DSZ ; decrement power. Done yet? GTO E ; no, keep going RTN ; yes, we're done

LBL E ; product so far is in X LSTx ; retrieve base (Y gets running product) MUL ; new product GTO D ; next iteration

Excuse the apparently arbitrary label selection. It nestles into a larger "suite". It handles floating point and integers (y**-x in integer mode returns zero with carry set. That seemed to be in keeping of the style of a 16C. Its major flaw is that it doesn't work for all real values of X--only integral ones. That sufficed for my purposes. I'm looking forward to seeing an implementation for all real X.

Also excuse my instruction mnemonics. They were etched on my brain long before I discovered that there was a "standard" notation system.

Cameron

      
Power notation syntax
Message #17 Posted by Andrés C. Rodríguez (Argentina) on 23 Oct 2001, 9:57 a.m.,
in response to message #1 by Cameron

Nice to see someone using y**x for power notation. At least for me, it has been a long time from seeing this (Fortran course in 1978), after that, y^x (Basic) became more usual... until modern word processors allow for superscript!


[ Return to Index | Top of Index ]

Go back to the main exhibit hall