WP34s: a small request

[fhub]

Would it be possible to change the new error message #20 from "Solve failed" to the more general message "no solution"?

Then it could also be used in many other programs (with ERR 20) whenever there is no solution for a problem. :-)

Franz

[Crawl]

Seems like a bad idea; that implies the solver knows for a fact that there is no solution, rather than that it merely failed to find a solution. "Routine failed" would be a more general message.

[fhub]

Ok, then let's forget it - I'll do it with an alpha message (which cost lots of bytes however :-(

[Walter B]

"Es war schon immer etwas teurer, einen besonderen Geschmack zu haben."

Old tobacco ad.

[fhub]

Quote:

---

"Es war schon immer etwas teurer, einen besonderen Geschmack zu haben."

Old tobacco ad.

---

Mein Gott Walter, /* dieses Lied kennst Du sicher, oder? */

Du bist heute wohl wieder besonders zum Scherzen aufgelegt! ;-)

[Paul Dale]

We also have a recently included three characters in two steps command. Not the easiest to use from the keyboard but it is there. The assembler will support it nicely however.

- Pauli

[fhub]

Yes I know, but I just thought that such a "no solution" message would be useful for many programs.

BTW, another question:

is it possible anyhow to switch on (in a program) the big "C" announciator as it is done automatically after a [CPX] operation?

I've tried it with 0 ENTER [cpx]+ at the end of a program before RTN (which doesn't change the results in the previous X and Y) to indicate a complex solution, but in a running program this doesn't seem to work: the "C" just didn't appear when the program stopped!?

I did it then with alpha commands, but that needs quite a few bytes ...

Franz

[Paul Dale]

You found the way: use alpha.

It doesn't use quite a few bytes BTW: three instructions I think. Clear alpha, a 'C' and one of the alpha views. The big 'C' is just a capital 'C'.

- Pauli

[fhub]

Quote:

---

three instructions I think. Clear alpha, a 'C' and one of the alpha views.

---

Yes, that's exactly as I did it. :-)

It was just for a little extension of the implemented SLVQ, which also return complex solutions:

/*

  Quadratic Equation:

  a*x^2 + b*x + c = 0

  Input:  a [ENTER] b [ENTER] c

  Output: X:x1 Y:x2, or X:Re(x0) Y:Im(x0) for complex solution (alpha:"C")

  Used:   Stack

*/

001 LBL'QE'

002 CL[alpha]

003 RCL/ Z

004 x[<->] Z

005 RCL+ X

006 /

007 +/-

008 x[^2]

009 RCL L

010 x[<->] Z

011 -

012 x<0?

013 SKIP 05

014 [sqrt]

015 RCL+ Y

016 x[<->]y

017 RCL- L

018 SKIP 05

019 ABS

020 [sqrt]

021 x[<->]y

022 [alpha] C

023 VW[alpha]+ X

024 RTN

Edited after Marcus' hint!

Edited: 13 Aug 2011, 4:50 a.m. after one or more responses were posted

[Marcus von Cube, Germany]

Franz, be aware that the following two lines are not the same:

[alpha]'C'

[alpha] C

The first line creates a multi byte string which is long enough to hold three characters in just two steps, the second one is just a single step. The effect on the alpha register is the same.

In order to enter a three letter alpha from the keyboard (only available in program mode) press f ENTER f ENTER. You'll see the single quote after [alpha] as a visual feedback.

[fhub]

Thanks Marcus - one more byte saved! :-)

If I optimize it even more (and get it down to 14 steps), then maybe Pauli will use it as code for SLVQ ... ;-)

Franz

[Paul Dale]

Figure out a compact way to implement this version and I'm not overly concerned about the number of steps :-)

The current SLVQ suffers from cancellation issues.

- Pauli

[fhub]

Quote:

---

Figure out a compact way to implement this version and I'm not overly concerned about the number of steps :-)

---

Well, although not much more code to do it this way, I don't think this would be worth the troubles (or extra bytes).

What then do with complex solutions (i.e. if b^2-4ac<0)? Make a complex division? Well, I guess this would destroy again the accuracy gained from these formulas.

And even with these special formulas you get completely wrong result for ill-conditioned equations - take for example:

x^2 + 10^8*x + 1 = 0

You'd need an internal accuracy of at least 16 digits to get correct solutions for this equation.

[Paul Dale]

Quote:

---

Well, although not much more code to do it this way, I don't think this would be worth the troubles (or extra bytes).

---

I think you're wrong :-) Numerical mathematics is a bizarre and arcane art quite unlike the rest of mathematics. This sort of attention to detail makes a difference.

Quote:

---

What then do with complex solutions (i.e. if b^2-4ac<0)? Make a complex division? Well, I guess this would destroy again the accuracy gained from these formulas.

---

I doubt the more complex formulas are losing accuracy here either :-) The accuracy is coming from the b - sqrt(...) where the sqrt(...) term is approximately b. i.e. b - b produces a catastrophic cancellation.

The b^2 - 4ac also risks cancellation but this is less critical. One does have to be careful about calculating this for b^2 approximately equal to 4ac however since it is easy to get the sign wrong :-( e.g. a simplistic application of a fused multiple add here can cause problem. Fortunately the 34S has a complex multiply that can achieve the desired calculation without any rounding issues.

Quote:

---

And even with these special formulas you get completely wrong result for ill-conditioned equations - take for example:

x^2 + 10^8*x + 1 = 0

You'd need an internal accuracy of at least 16 digits to get correct solutions for this equation.

---

Clearly, this is exceeding the precision available. But there are plenty of cases in the middle where b^2 >> 4ac and the naïve formula results in b - almost(b) with a huge loss of precision and a resulting incorrect answer.

This middle region is where these more complicated formulas come into their own. Try this reference for a better description.

- Pauli

Edited: 13 Aug 2011, 6:22 a.m. after one or more responses were posted

[fhub]

Quote:

---

Numerical mathematics is a bizarre and arcane art quite unlike the rest of mathematics.

---

Correct, and that's exactly the reason why I don't like it very much. ;-)

And if I really need such an excessive accuracy - well, them I'm using a CAS with almost _arbitrary_ precision. :-)

[Paul Dale]

The 34S needs to produce correct results and retaining accuracy is mandatory. Look at the trouble I had with the statistical functions. With Dieter's assistance, we've managed to get them accurate with correct rounding.

As I mentioned in the discussion about the 35s, a calculator must produce correct results all the time. It is a tool and the tool must work properly.

- Pauli

_{edit to fix silly grammatical / typo error -- I'm having a bad evening :-(}

Edited: 13 Aug 2011, 6:24 a.m. after one or more responses were posted

[fhub]

Quote:

---

The 34S needs to produce correct results and retaining accuracy is mandatory.

---

Ok, I'll work on these 2 (more accurate) formulas! ;-)

But getting also the complex case will certainly need a few bytes more ...

[Paul Dale]

A quadratic solver using the stable formulas which also handles complex roots is pretty much guaranteed to be included. I rate correct and stable answers as more important than everything else :-)

The reason the current solver is there is because it is very small, neat and concise. It is a lovely example of keystroke programming. That it suffers from instability is a minus that should be corrected in due course (sorry Gerson).

- Pauli

Edited: 13 Aug 2011, 6:21 a.m.

[fhub]

Quote:

---

A quadratic solver using the stable formulas ...

---

Well, at least for real solutions I already have the code - it's 21 steps compared with 14 in xrom.c.

But the complex case needs more space (and time ... ;-))

[Marcus von Cube, Germany]

Quote:

---

But the complex case needs more space (and time ... ;-))

---

In the complex case the cancellation isn't a problem except for the term below the square root. The subtraction b - sqrt(..) is never carried out. :-)

[fhub]

Quote:

---

In the complex case the cancellation isn't a problem except for the term below the square root. The subtraction b - sqrt(..) is never carried out. :-)

---

Correct - and right I'm finished with it! :-)

With LBL and RTN it's 37 steps (compared to 16 in xrom.c), could be made 3-4 steps shorter by removing 2 tests (like e.g. a=0?) but this would give 'strange' error messages for some cases.

Two things are not optimal:

- I needed an extra register to get it working also in SSIZE4 mode (I've used "I")

- after execution the solution is in X and Y, but Z and T contain 'garbage' (intermediate results): avoiding this would need quite a few more bytes or at least other Rnn registers, and I don't think this would be a good idea

So let me know if these circumstances are ok and I should post my solution, or if I can throw it away ... ;-)

Franz

[Walter B]

Quote:

---

Two things are not optimal:

- I needed an extra register to get it working also in SSIZE4 mode (I've used "I")

---

That's suboptimal indeed, since "I" is the second half of complex last x. How about taking "A" instead? Another opportunity may be confining this function to stack size 8 (I knew we would need this stack size for certain tasks d:-)

Quote:

---

- after execution the solution is in X and Y, but Z and T contain 'garbage' (intermediate results)

---

Applying your measure of tidyness (of an earlier thread), the stack levels not needed for the result shall contain what was on the stack before calling complex SLVQ.

Just my 20m€,
Walter

[fhub]

Quote:

---

That's suboptimal indeed, since "I" is the second half of complex last x.

---

And where's the problem???

"I" is not used at all in this calculation, and also "L" contains garbage at the end - what counts are the 2 solutions in X and Y, nothing else!

Quote:

---

Applying your measure of tidyness (of an earlier thread), the stack levels not needed for the result shall contain what was on the stack before calling complex SLVQ.

---

How could I retain 3 coefficients (a,b,c) AND return 2 solutions (x1,x2) on a 4-level stack, can you please explain this?

So the best solution would be _replacing_ the 3 coeffs by the 2 solutions - but this would need again (not so few) bytes.

Ok, I've now finally lost my interest in trying to fulfill your extremely high 'requirements' - if you want and need all these conditions then do it yourself. :-(

Franz

[Walter B]

Quote:

---

"I" is not used at all in this calculation, ...

---

versus

Quote:

---

I needed an extra register to get it working also in SSIZE4 mode (I've used "I")

---

Now what??

Quote:

---

... what counts are the 2 solutions in X and Y, nothing else!

---

I don't get it. How do you place 2 complex solutions in two stack levels :-?

[fhub]

Sorry, don't you really get it or do you just want to annoy me?

"I" is NOT used for CALCULATION but just for STORAGE!

Quote:

---

I don't get it. How do you place 2 complex solutions in two stack levels :-?

---

Since these 2 solutions are complex conjugates, it's enough to put Re(x0) in X and Im(x0) in Y.

The 2nd solution is just Re(x0)-Im(x0)*i

Edited: 13 Aug 2011, 2:09 p.m.

[Walter B]

Thanks for your explanation of the 2 solutions. I should have known that but I didn't :-/

For the other point:

Quote:

---

"I" is NOT used for CALCULATION but just for STORAGE!

---

You're shouting this - but you didn't tell this before. Keep cool, boy 8-) "Needing" a register may be for any reason (in a calculator I guess most probably for calculations but - as stated several times already - my abilities in clairvoyance are limited). No, I didn't want to annoy you, though I admit it's very easy to tease you - and easy to foresee you will declare your resignation (once again) but won't do it ;-)

[fhub]

Quote:

---

You're shouting this

---

No, I didn't shout - it's just faster to capitalize words that I want to mark than putting them between 2 "bold"s.

Quote:

---

No, I didn't want to annoy you, though I admit it's very easy to tease you

---

Ahhh, interesting, so you admit you want or try to tease me ... ;-)

Edited: 13 Aug 2011, 6:43 p.m.

[Paul Dale]

Quote:

---

How about taking "A" instead?

---

No no no no no!!!!!!! Lettered things are special and shouldn't be used in libraries.

Fortunately, this is intended for XROM. With a suitable prologue and epilogue added to the code, there are another five registers available. These are numbers 00 through 04 and they are completely separate from the same numbered registers the user sees. There are also fifteen extra flags 00-14 FWIW. XROM also operates with a stack size of four regardless of what the user sets -- this might cause problems with stack drop.

So, what I'd be doing is to preserve T & I in two of these extra registers, use a third for the temporary storage of the value. A fourth preserving X if you want Last X done like an internal operation.

Then do the calculation and finally fix things up at the end. Yes this will take more steps. Often quite a few more.

For the b^2 - 4ac calculation use a complex multiply.

Assuming a, b & c are in X, Y & Z respectively (also should work for c, b & a in X, Y & Z).

    4          a b c T

    *          4a b c T

    RCL Y      b 4a b c

    cpx *      b^2-4ac b(4a+c) b c

    x<>y       b(4a+c) b^2-4ac b c   if only we'd kept the DROPY instruction :-)

    DROP       b^2-4ac b c c

Which, assuming I've got that correct, will calculate the difference without any loss of precision. Also, 4a is still in I.

- Pauli

Edited: 13 Aug 2011, 8:32 p.m.

[Paul Dale]

I noticed your code multiplies 2a by 2c which makes more sense. The principle stays the same luckily.

- Pauli

Edited: 13 Aug 2011, 9:23 p.m.

[Marcus von Cube, Germany]

The register isn't a problem if your code makes it to XROM because we have some private registers there for exactly this purpose. I wouldn't mind that the stack gets used, other functions do the same. It's just not a standard function but a routine. With the help of the private registers we may even save stack level T and restore it. I have no idea what would be the "correct" result for L in this case.

Please post your solution! Chances are high that it will replace what we have now.

[fhub]

Well, with additional internal registers many thing would of course be much easier to code. And about "L" I won't care at all - after such a lot of operations there's just no useful value which should intentionally be placed in "L".

I've now checked all 8 test equations at the link of Gerson, and for 7 of them I get exactly the same solutions. Only in Case 5 my program returns 2 equal real solutions (the same as the real part on this site), when in fact there is a imaginary part, but it is about 3e-10. For this special equation (Case 5) both expressions b^2 and 4ac are VERY big and differ only at the 17th or 18th digit, so that's of course far beyond the 34s precision.

Although I would have an idea how to solve even this problem (b^2-4ac could be factorized as (b+2*sqrt(ac))*(b-2*sqrt(ac)) which avoids those big numbers and would certainly work), but I don't think it would be worth the troubles and many extra bytes (again an extra codeblock if ac<0!). But also this method could not guarantee always correct solutions if the numbers are still bigger.

So here's my code:

/*

  Quadratic Equation (extended precision)

  a*x^2 + b*x + c = 0

  Input:  a [ENTER] b [ENTER] c

  Output: X:x1 Y:x2, or X:Re(x0) Y:Im(x0) for complex solution (alpha:"C")

  Used:   Stack, I

*/

001 LBL'QE'

002 RCL+ X

003 x[<->] Z

004 x=0?

005 ERR 18

006 RCL+ X

007 x[<->]y

008 STO I

009 x[^2]

010 RCL Y

011 RCL[times] T

012 -

013 x<0?

014 SKIP 13

015 [sqrt]

016 RCL I

017 x<0?

018 SKIP 03

019 x[<->]y

020 +/-

021 x[<->]y

022 -

023 x[!=]0?

024 STO/ Z

025 x[<->]y

026 /

027 SKIP 09

028 ABS

029 [sqrt]

030 RCL/ Y

031 RCL I

032 RCL/ Z

033 +/-

034 CL[alpha]

035 [alpha] C

036 VW[alpha]+ X

037 RTN

[Paul Dale]

Steps 19, 20 and 21 are a complex conjugate, so replacing them with CONJ saves two more steps. Then the conditional can be reversed and the skip removed for a third step saved.

That is replace:

    017 x<0?

    018 SKIP 03

    019 x[<->]y

    020 +/-

    021 x[<->]y

with

    017 x>0?

    018 CONJ

and fix the skips up.

- Pauli

[fhub]

Quote:

---

Steps 19, 20 and 21 are a complex conjugate, so replacing them with CONJ saves two more steps.

---

A nice trick indeed! (maybe you should make a doc file "tips & tricks for 34s programming" ;-))

But this CONJ also changes RegI (where I have stored the coefficient b), so I hope you haven't overlooked this.

Quote:

---

there are another five registers available. These are numbers 00 through 04 and they are completely separate from the same numbered registers the user sees

---

Hey, that's quite unfair if you have more registers than we have! ;-)

Why can't we use them in our programs too? That would simplify programming a lot ... :-)

Franz

[Paul Dale]

Quote:

---

Hey, that's quite unfair if you have more registers than we have! ;-)

Why can't we use them in our programs too? That would simplify programming a lot ... :-)

---

You can have them if you're willing to forgo solve, integrate, sum, product and the derivatives. Oh and now the quadratic solver. I'm sure most people would like to keep these functions and not have the extra registers.

It is also not a case of having an extra five registers. All the above are limited to using no more than these five registers in total. They don't even get the luxury of having the stack preserved across callbacks to user code.

- Pauli

[Paul Dale]

Quote:

---

A nice trick indeed! (maybe you should make a doc file "tips & tricks for 34s programming" ;-))

---

I'd like to eventually. I've got other things for the moment.

Quote:

---

But this CONJ also changes RegI (where I have stored the coefficient b), so I hope you haven't overlooked this.

---

I stopped using register I and restore it afterwards :-)
I'm using a complex multiply to retain accuracy under the square root which also loses I.

- Pauli

[fhub]

Quote:

---

I stopped using register I and restore it afterwards :-)
I'm using a complex multiply to retain accuracy under the square root which also loses I.

---

Ok, then it's no problem.

BTW, since you mention accuracy:

What do you think about using (b+2*sqrt(ac))*(b-2*sqrt(ac)) instead of b^2-4ac ?

This would eventually help when b^2 and 4ac are VERY big and their difference is beyond the 34s precision - I guess this could also solve example 5 in Gerson's link!?

OTOH it would require quite a few bytes more, so I'm not sure if it's really worth to be used.

And one more thing (rather for Walter than for you):

If I have understood your code in xrom.c, then SLVQ skips the next step (in a program) for complex solutions but NOT for real solutions, right?

This should be described clearly in the manual, because it's not so obvious.

Franz

Edited: 14 Aug 2011, 8:04 a.m.

[Paul Dale]

Quote:

---

What do you think about using (b+2*sqrt(ac))*(b-2*sqrt(ac)) instead of b^2-4ac ?

---

It should address the overflow issues, or at least push them off a bit. I'd be concerned about the second term when b is approximately 2 sqrt(ac). It might not be an issue, but it would have to be investigated.

Quote:

---

This would eventually help when b^2 and 4ac are VERY big and their difference is beyond the 34s precision - I guess this could also solve example 5 in Gerson's link!?

---

The current code handles example 5. The complex multiply addresses the b^2 - 4ac precision loss concerns. The internal code evaluates the multiplies exactly and does the subtraction without any rounding or truncation until the result is returned. i.e. it should be as good as is possible without resorting to extended precision.

- Pauli

[fhub]

Quote:

---

I'd be concerned about the second term when b is approximately 2 sqrt(ac). It might not be an issue, but it would have to be investigated.

---

Well, but anywhere you would have to stop - what should you do when those 2 terms are almost equal? You can't completely avoid subtractions in any cases ...

Quote:

---

The current code handles example 5. The complex multiply addresses the b^2 - 4ac precision loss concerns.

---

Aaah, that's fine - so this problem seems to be solved at least with this extended precision. Of course we can't never handle ALL situations correctly when the numbers are getting to big.

Franz

[Paul Dale]

Much goodness in subversion. No build yet.

- Pauli

[Gerson W. Barbosa]

Quote:

---

A quadratic solver using the stable formulas which also handles complex roots is pretty much guaranteed to be included.

---

That's really good news!

Palmer O. Hanson has written an excellent quadratic solver that handles some really tough examples:

A "Cadillac " Quadratic Solver for the hp 33s and HP35s

Gerson.

[Walter B]

Quote:

---

Quote:

---

Numerical mathematics is a bizarre and arcane art quite unlike the rest of mathematics.

---

Correct, and that's exactly the reason why I don't like it very much.

---

Bad luck - calculators are numerical mathematics all over ;-)

[fhub]

Quote:

---

Bad luck - calculators are numerical mathematics all over ;-)

---

That's why I'm USING but not PRODUCING calculators! ;-)

[Gerson W. Barbosa]

Quote:

---

Figure out a compact way to implement this version

---

Here's my first attempt at it. Like Franz's, it handles Palmer's examples (except number 5). Too many steps though, "cut a few and it'll be perfect" :-)

Gerson.

// Source file(s): wp34s-0.dat

001 LBL B

002 CL[alpha]

003 x[<->]y

004 x[>=]0?

005 SKIP 06

006 ENTER[^]

007 SIGN

008 STO[times] Y

009 STO[times] Z

010 STO[times] T

011 R[v]

012 RCL Z

013 x=0?

014 ERR 18

015 +/-

016 STO+ X

017 STO 00

018 STO+ X

019 RCL[times] Z

020 RCL Y

021 x[^2]

022 +

023 x<0?

024 SKIP 08

025 [sqrt]

026 +

027 RCL/ 00

028 R[^]

029 RCL/ L

030 STO+ X

031 +/-

032 SKIP 07

033 +/-

034 [sqrt]

035 RCL/ 00

036 +/-

037 x[<->]y

038 RCL/ 00

039 [alpha] C

040 VW[alpha]+ X

041 RTN

// 41 total instructions used.

// 41 total words used.

// 41 single word instructions.

// 0 double word instructions.

_{Edit to correct the listing as per Franz's observation below.}

Edited: 14 Aug 2011, 6:31 p.m. after one or more responses were posted

[fhub]

Quote:

---

006 ENTER[^]

007 SIGN

008 STO[times] X

009 STO[times] Y

010 STO[times] Z

011 R[v]

---

I guess these 3 "STO[times]" should be Y,Z and T instead.
With X,Y and Z in fact none of these values is changed.

[Gerson W. Barbosa]

You're right, thanks! It's correct on my wp34s, but not on the emulator. Will fix the listing.

Gerson.

[Paul Dale]

The code I've put into xrom is:

81ee: 5b59       LBL 89         Entry: QUADRATIC SOLVER

81ef: 5d28       XEQ 40         -> Common entry code

81f0: 3f00       STOS 00

81f1: 2b6d       RCL I

81f2: 2404       STO 04

81f3: 0109       R[v]

81f4: 2c64       RCL+ X

81f5: 3266       x[<->] Z

81f6: 0017       x=0?

81f7: 2312       ERR 18

81f8: 2c64       RCL+ X

81f9: 2402       STO 02

81fa: 2b65       RCL Y

81fb: 0603       [cmplx][times]

81fc: 2401       STO 01

81fd: 001a       x<0?

81fe: 530b       SKIP 11

81ff: 020a       [sqrt]

8200: 0107       x[<->]y

8201: 0109       R[v]

8202: 001c       x>0?

8203: 0536       CONJ

8204: 0302       -

8205: 0018       x[!=]0?

8206: 2865       STO/ Y

8207: 2f02       RCL/ 02

8208: 5d01       XEQ 01

8209: 5e29       GTO 41         -> Normal exit code

820a: 0205       ABS

820b: 020a       [sqrt]

820c: 2f02       RCL/ 02

820d: 2b66       RCL Z

820e: 2f02       RCL/ 02

820f: 0003       +/-

8210: 5d01       XEQ 01

8211: 019c       TOP?

8212: 5301       SKIP 01

8213: 5e2a       GTO 42

8214: 0171       CL[alpha]

8215: 4143       [alpha] C

8216: 9d64       VW[alpha]+ X

8217: 5e29       GTO 41         -> Normal exit code

8218: 5b01       LBL 01

8219: 2b00       RCL 00

821a: 246c       STO L

821b: 2b04       RCL 04

821c: 246d       STO I

821d: 0108       [cmplx]x[<->]y

821e: 2b04       RCL 04

821f: 2b01       RCL 01

8220: 0108       [cmplx]x[<->]y

8221: 013a       RTN

Label 40 is the standard entry code -- checks for nesting in xrom.

Label 41 is basically a RTN.

Label 42 is a RTN+1.

It deals with all of Palmer's examples properly.

- Pauli