More trig--specifically speed vs accuracy



#18

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?


#19

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.


#20

Bruce,

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


Viktor


#21

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


#22

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


#23

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


#24

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.

#25

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 ...


#26

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.


#27

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.

#28

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


#29

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.


#30

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


#31

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.


#32

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.


#33

...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

#34

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!


Possibly Related Threads...
Thread Author Replies Views Last Post
  48G vs 49G+ User RPL Speed Comparison John Colvin 7 321 11-16-2013, 10:07 PM
Last Post: Han
  Trig vs hyperbolic handling differences in Prime CAS Michael de Estrada 3 239 11-08-2013, 06:26 PM
Last Post: Mark Hardman
  Trig Functions Howard Owen 11 541 09-16-2013, 02:53 PM
Last Post: Fred Lusk
  trig scales on the Post Versalog slide rule Al 12 566 09-15-2013, 06:01 AM
Last Post: John I.
  WP-34S: Speed of y^x Marcel Samek 1 175 09-14-2013, 07:31 PM
Last Post: Paul Dale
  WP-34S function execution speed ? Gene Wright 4 253 09-04-2013, 05:40 PM
Last Post: Paul Dale
  How much accuracy does one actually need? Matt Agajanian 23 784 08-26-2013, 12:46 PM
Last Post: Kimberly Thompson
  Some trig help, please Matt Agajanian 33 1,028 04-12-2013, 06:39 PM
Last Post: Matt Agajanian
  Trig Identity question Namir 13 511 03-04-2013, 07:45 PM
Last Post: Eric Smith
  HP-39gII speed Mic 2 216 02-24-2013, 05:55 PM
Last Post: Thomas Klemm

Forum Jump: