HP Forums

Hi everybody, here we go (again):

Say you've just got that powerful, new HP calculator for your collection (say an HP33S
or some HP48/49 model) and you're itching to test it on some interesting problem.
Here's one of my own devising for you to try. Apart from being interesting in its own right (IMHO), it might actually be useful ! Consider for instance its use by someone teaching math to prepare interesting problems for his/her students (i.e.: K = 7.041776 anyone ? :-) Read on ...

Before we begin ...

We all know Pi is a transcendental number. By definition, this means Pi cannot be the root
of any polynomial equation with integer coefficients, whatever its degree (this is also
true of all other transcendental numbers such as e, sin(1), log(2), etc). In particular,
Pi cannot be a root of a cubic equation with integer coefficients, such as

        A*x^3 + B*x^2 + C*x + D = 0

The challenge:

"Write a program for your preferred HP calculator to find the integer coefficients A,B,C,D of the cubic equation:

        A*x^3 + B*x^2 + C*x + D = 0

which has a root that best approximates a given positive constant K (say, K = Pi), where the
coefficients are less or equal in absolute value to a maximum given integer value, M
(say, M = 9)."

Notes:

Without loss of generality, the A coefficient will always be positive
(i.e: it goes from 1 to M, both included) while the B,C, and D coefficients can be
either positive or negative (i.e: they go from -M to +M, both included). Your
program must input K (say Pi), M (say 9) and a very loose initial interval for the
root (say [3.1, 3.2]) and must output the coefficients A,B,C,D of the "best"
equation, the value X of the root, and the absolute difference between the
root X and the given constant K.

For instance, given the inputs K=Pi, M=20, Interval= [3.1, 3.2] it might
output something like (not the best solution, of course):

       A= 1, B= -10, C= 19, D= 8, X= 3.14162732179, DIF= 0.00003466820 

You are expected to deliver a general program, that will work for any positive constant K, maximum coefficient M and given initial interval [X1,X2], but for the purposes of testing the efficiency and accuracy of your solution you should try these two cases:

• Case 1: constant K = Pi, Max. coefficient M = 9, initial interval [3.1, 3.2]

• Case 2: same but with Max. coefficient M = 20
  

Hints:

• Absolutely refrain from using a PC. You'll learn nothing about writing and optimizing
  programs for your HP calculator and will ruin the pleasure of seeing your efforts
  and ingenuity succeed with your little HP machine. This is intended for real programmers using real HP calcs, not the "brute-force-using-a-zillion-Ghz-PC" kind.

• Though most any HP model can be used to successfully solve this challenge, it is advisable
  that you use a suitably fast HP calculator, say any model using a Saturn CPU (HP33S, HP32S,
  HP32SII, HP42S, HP48/49, HP-71B, etc), in order to attain reasonable running times.

• Even so, you'll need to exercise your ingenuity to achieve bearable running times.
  A brute-force approach is doomed to failure. Consider this:

  - for M = 9, your program will have to choose among no less than 61,731 possible cubic equations, from (1,-9,-9-,9) to (9,9,9,9)

  - for M=20, your program will have to choose among more than 1 million possible equations,
  (namely 1,378,420) from (1,-20,-20,-20) to (20,20,20,20)

• For the case constant K = Pi, Max. coefficient M = 9, initial interval [3.1, 3.2],
  it takes 7 min. 16 seconds to find the solution among the 61,731 cubic equations
  possible, output it and stop.

  
  

• For the case constant K = Pi, Max. coefficient M = 20, initial interval [3.1, 3.2],
  it takes 3 hour 4 min. to find the solution among the 1,378,420 cubic equations
  possible, output it and stop.

Now's your turn. Just put that brand new shiny HP33S, that wonderful classical HP42S, that ultrapowerful HP49, or that venerable HP-71B to the task, plus your ingenuity and RPN/RPL/BASIC/FORTH/Assembler programming skills of course. After all, if the good old HP-71B can do it in 3 hours using BASIC,
what incredible times will you be able to achieve with much faster machines using RPN/RPL, say ? Emulators very
welcome too, of course.

Best regards from V.

I've written a program for my HP

a real one (HP-32SII)

I found a result for M=8 as well as for M=20. The latter took about an hour for each value of A, evaluating B,C&D.
(think I need new batteries now...)

I must confess however, that I didn't fulfill the assignment to the letter; I haven't "rooted" the equation for each set of parameters. I know I make a principle error here, but I think the final results will turn out the same.



Am I supposed to put my answers here?

Hi Bram:

Bram posted:

"I've written a program for my HP a real one (HP-32SII) I found a result for M=8 as well as for M=20. [...] I think the final results will turn out the same.
Am I supposed to put my answers here? "

Yes, the results *and* your program for the 32SII, if you don't mind. That way we'll be able to see if your computed best solutions agree with mine and fully appreciate the programming ingenuity you used for the task, I'm pretty sure many people will be delighted to have a look at it.

As for the times, I understand it took your program some 20+ hours for the M=20 case, right ?

Thanks for your interest and best regards from V.

Ok, here's my go at it.

It's an HP71 program, 9 lines, 344 bytes.

The print statement alone takes up 72 bytes.

I didn't use an input statement to input K and M, but just inserted them in the program

I also left in some statements to compute the time and error as it goes along.

It could have been made shorter by leaving those out.

For K=PI and M=9, it takes 2 min, 44 sec to find A=1, B=-1, C=-8, D=4, err=-.00025653792

For K=PI and M=20, it takes 26 min, 52 secs to find A=6, B=-16, C=-15, D=19, err=-.00000160749



10 K=PI @ M=9 @ E=MAXREAL @ T9=TIME

20 FOR P=1 TO M @ T0=P*K

30 FOR Q=-M TO M @ T1=(T0+Q)*K

40 FOR R=-M TO M @ T2=(T1+R)*K

50 IF ABS(T2)>M+.5 THEN 70

60 T=ABS(FLOOR(T2+.5)-T2) @ IF T<E THEN E=T @ A=P @ B=Q @ C=R @ D=FLOOR(.5-T2) @ PRINT E

70 NEXT R @ NEXT Q @ NEXT P @ S=(3*A*K+2*B)*K+C @ X=K @ Z=MAXREAL

80 Y=((A*X+B)*X+C)*X+D @ IF Y=Z THEN 90 ELSE Z=Y @ X=X-Y/S @ GOTO 80

90 PRINT "A=";A;"B=";B;"C=";C;"D=";D;"X=";X;"DIF=";X-K;"TIME=";TIME-T9;"SECS"



The program even works for M=30, but the simple little root finder in line 80

gets stuck in an infinite loop because the computed value of Y is dominated by round-off error

Replacing 'IF Y=Z THEN 90' with 'IF Y<1E-9 THEN 90' gets out of the infinite loop

but returns a less accurate root

Hi, Rodger:

Rodger posted:

"Ok, here's my go at it.
It's an HP71 program, 9 lines, 344 bytes. [...]For K=PI and M=9, it takes 2 min, 44 sec to find A=1, B=-1, C=-8, D=4, err=-.00025653792 --- For K=PI and M=20, it takes 26 min, 52 secs to find A=6, B=-16, C=-15, D=19, err=-.00000160749"

Congratulations ! Your solutions are absolutely correct and your program is as short and fast as it gets, very good work indeed !

Just for the record, the solutions are:

  M = 9

     x3 - x2 - 8x + 4 = 0

        x = 3.14133611565 (0.00025653794)

  M = 20

     6x3 - 16x2 - 15x + 19 = 0

        x = 3.14159104610  (0.00000160749)

Just in case it might interest you, here is an assortment of results for larger values of M:

  M = 49

     19x3 - 47x2 - 30x - 31 = 0

        x = 3.14159235658  (2.97E-07)

  M = 99

     33x3 - 92x2 - 65x + 89 = 0

        x = 3.14159264441 (9.18E-09)

  M = 199

     37x3 - 114x2 - 36x + 91 = 0

        x = 3.14159265383  (2.39E-10 )

Best regards from V.

By now I know that my answers are right and here's how I got them.

I already announced that I kind of cheated. Considering that the root finding part will take a lot of computing time, I decided to skip it by evaluation the function for point K for every combination of ABCD within the given limits. The combination for which the function is closest to zero, the real root will most probably be closest to K. "Most probably", because it isn't necessarily so, but as we have integer coefficients and hence the function is kind of stepping through its possible values, I thought it very likely to get the right answers all the same.

During the process every ABCD evalution that is better than the best so far, is set apart.

I realized this curious thinking in my program

When you follow the computing progress, there are times that you would like to say to the machine: "hey, don't evaluate this part any further. It's no use". Those are the parts where the program could be speeded up. F.e. when for any combination of ABC the function is already to much positive for D=-20, then you may skip the entire D-range. (I may squeeze this speed-up into the program some day).

When I finally found an ABCD combination, I did compute the real root of the function using the Solve program of my HP-32SII (EQN: D+X*(C+X*(B+X*A)) )(I hate brackets ;-) )

Then I subtracted Pi

and it appeared I got the right answers:

1,-1,-8,4 root=3.14133611565 err= -.000252653794

6,-16,-15,19 root=3.14159104610 err=-.00000160749



Indeed, the latter took almost 24 hours to complete, but it wwas a vary nice day.

Hi Bram,

Bram posted:

"By now I know that my answers are right and here's how I got them. [...] The combination for which the function is closest to zero, the real root will most probably be closest to K. "Most probably", because it isn't necessarily so, but as we have integer
coefficients and hence the function is kind of stepping through its possible values, I thought it very likely to get the right answers all the same."

You're right, very good idea. This kind of shows how plausible heuristics can go a long way toward making feasible and successful an otherwise hopeless brute-search approach.

"... when for any combination of ABC the function is
already to much positive for D=-20, then you may skip the entire D-range."

Right, too. Another very good heuristic. This speeds the program immensely by essentially removing the innermost loop and thus halving (or better) processing time.

"Then I subtracted Pi and it appeared I got the right answers:"

You certainly did.

"Indeed, the latter took almost 24 hours to complete, but it was a vary nice day. "

I'm glad you took some enjoyment from my 'challenge', thanks for submitting your clever program to this forum, for using up so much juice from your batteries, and last but not least, for your interest in my humble ramblings.

Best regards from V.

Valentin,

Have you enough results (or can you generate them easily) to indicate whether the results get monotonically better (i.e. more and more accurate) with larger and larger values of M?

Actually, the answer must be NO, since for M=199, your biggest coefficient is only 114, so this must also have been the best solution for M = 115 . Similarly, for your other large-M results, the answers do not bump up to the limit set by M.

To what extent are these solutions local minima? i.e. do some of these solutions repeat for other values of M? I would think so, again based on your M = 199 solution. If I understand the challenge correctly, this solution must be the best FOR ALL VALUES OF M FROM 115 to 199. Otherwise, your answer would have a larger coefficient than 114 in it.

How far beyond 199 do you have to go until there is the next better solution?!?!

Hi, Dave:

Dave posted:

"Actually, the answer must be NO, since for M=199, your biggest coefficient is only 114, so this must also have been the best solution for M = 115."

That's correct ...

"To what extent are these solutions local minima? i.e. do some of these solutions repeat for other values of M? I would think so, again based
on your M = 199 solution. If I understand the challenge correctly, this solution must be the best FOR ALL VALUES OF M FROM 115 to
199. Otherwise, your answer would have a larger coefficient than 114 in it."

Again correct. But please, don't shout ! :-)

"How far beyond 199 do you have to go until there is the next better solution?!?!"

I didn't test it extensively, but the solution for M=114 is good until M=255 (inclusive), at least. Also, a run with M=360 gives

   ---------------------------------------------------------

      M    a    b    c    d            x      Diff=ABS(Pi-x)

    ---------------------------------------------------------

    360  132 -319 -229 -225    3.14159265348  0.00000000011

       132x3 - 319x2 - 229x - 225 = 0

        x =  3.1415926534831+

    For K = e = 2.71828182846+

    ---------------------------------------------------------

      M    a    b    c    d            x      Diff=ABS(e-x)

    ---------------------------------------------------------

      9    3   -5   -6   -7    2.71823262630  0.00004920216

     49    6   -4  -25  -23    2.71828239152  0.00000056306

     99   29  -63  -64   57    2.71828183050  0.00000000204 

    255   75 -179  -15 -143    2.71828182831  0.00000000015 

The last line features a solution valid from M = 179 to at least 255, and the resulting equation:

       75x3 - 179x2 - 15x - 143 = 0

        x =  2.71828182831+

       4x3 - 26x2 - 11x - 30 = 0

        x = 7,041776+ . 

Does it ring some bells ? :-) It would make a very fine, patriotic assignment for any US students to solve, right ? It could also be used to construct some puzzle-like problem where after telling some story a number of conditions are given resulting in this equation, then its real root is precisely that most famous date.

Thanks for your thoughtful insights and your interest, and

Best regards from V.

Edited: 28 May 2004, 5:26 a.m.

Hi,

Just some optimizations on the same topic:

With UTIL/1 binary (AWRITE) and an HP86B:

10 K=PI @ M=9 @ N=M+1 @ O=N*K @ P=N/K @ E=1.E99 @ CLEAR

20 FOR A=1 TO M @ T0=A*K @ IF ABS(T0)*K*K-N>O*(K+1) THEN 100

22 AWRITE 0,0 @ DISP A; T0

30 B0=MAX(-M, INT((-P-M)/K-T0+0.5)) @ B1=MIN(M, INT(P+M)/K-T0+0.5)) @ IF B0>B1 THEN 100

32 AWRITE 0,40 @ DISP B0; B1

40 FOR B=B0 to B1 @ T1=(T0+B)*K @ IF ABS(T1)*K-N>O THEN 90

42 AWRITE 1,0 @ DISP B; T1

50 C0=MAX(-M, INT(-P-T1+0.5)) @ C1=MIN(M, INT(P-T1+0.5)) @ IF C0>C1 THEN 90

52 AWRITE 1,40 @ DISP C0; C1

60 FOR C=C0 TO C1 @ T2=(T1+C)*K @ D=-INT(T2+0.5) @ IF ABS(D)>M THEN 80

62 AWRITE 2,0 @ DISP C; T2

70 IF ABS(T2+D)<E THEN F=A @ G=B @ H=C @ I=D @ E=ABS(T2+D) @ AWRITE 5,0 @ DISP F; G; H; I; E

80 NEXT C

90 NEXT B

100 NEXT A @ AWRITE 6,0 @ END

First it avoid looping on non-feasible A, B and C parameters.

Second it doesn't really solve the polynom, it can perhaps find a minimum.

Otherwise it's quite efficient (you can remove lines 22, 32, 42, 52, 62).

Bye,

Olivier

Edited: 29 May 2004, 6:57 a.m.

I've speeded up the program. I thought it wouldm't make a big difference, but it surely does. Now I get (the same) answers in only 3 hours in stead of almost 24.
The new program is with 3 lines more code significantly faster. You'll find the extension in the D-labeled part.

Thanks for your interest and program, Olivier, but not having an HP-86 at hand (and I suspect most MoHP regulars don't either) it would be nice if you could provide sample results and timing, to see that it delivers the correct results and also how it does compare with previous solutions.

Best regards from V.

Sorry, Bram, but the link to your new version doesn't work.

Best regards from V.

Would you be so kind as to mail me?
Drop the obvious in my account: pitquim@yahoo.es.quita.esto

Thanks in advance

Raul L

Edited: 31 May 2004, 10:02 a.m.

Best regards from V.

I am very very much ashamed.
Please retry

Gracias

Hi,

Sorry for the delay:

For the HP86B it take about 50 sec for M=20, and 29700sec for M=200 with a better optimized program.

Anyway here is a version for an HP49G+ with all the optimizations:

%%HP: T(3)A(R)F(.);

\<< \pi SWAP MAXR OVER 1. + DUP 5. PICK DUP2 * UNROT / DUP 6. PICK + 7. PICK / { } \-> K M E N O P Q F

  \<< CLLCD K 1. M

    FOR A A 1. DISP

      IF DUP ABS K SQ * N - O K 1. + * \<=

      THEN DUP M NEG OVER Q + NEG 0. RND MAX DUP UNROT + K * SWAP M Q 4. PICK - 0. RND MIN

        IF DUP2 \<=

        THEN

          FOR B

            IF DUP ABS K * N - O \<=

            THEN DUP M NEG OVER P + NEG 0. RND MAX DUP UNROT + K * SWAP M P 4. PICK - 0. RND MIN

              IF DUP2 \<=

              THEN

                FOR C DUP 0. RND

                  IF DUP ABS M \<=

                  THEN

                    IF DUP2 - ABS E <

                    THEN DUP2 - ABS DUP 6. DISP 'E' STO DUP NEG C B A 4. \->LIST REVLIST DUP 'F' STO 7. DISP

                    END

                  END DROP K +

                NEXT

              ELSE DROP2

              END DROP

            END K +

          NEXT

        ELSE DROP2

        END DROP

      END K +

    NEXT DROP E F OBJ\-> \->ARRY DUP 'S' STO

    \<< S X PEVAL

    \>> 'X' K ROOT K

  \>>

\>>

It take 5sec for M=9 and 40sec for M=20. Sorry I didn't test other M.

Olivier

P.S. all the results are ok with yours.

Edited: 1 June 2004, 2:08 a.m.

No luck. This time it says "Service Unavailable".

May I suggest you simply include the whole listing in a normal post, between opening and closing "pre" codes
to preserve the formatting ? It would be much simpler ...

Thanks for your efforts and

Best regards from V.

Thank you very much, Olivier ! Magnificent effort !!

I see you master a jaw-dropping range of machines, from such an old timer as an HP-86 to the very latest HP49G+. Amazing indeed !

I did write a lot of commercial engineering programs for the HP86 and HP87 machines back in their time. Most of them were written in their built-in dialect of BASIC, but laced with a number of binary programs I wrote myself, thanks to the Assembler ROM, lots of documentation and the wonderful fact that you could have as many as 5 binary programs loaded at a time versus just one in the former HP85 models.

It was great fun when your new keywords did work ok. I specially remember a big binary I wrote which accomplished all kinds of matrix operations (at assembler speeds) but with very large, sparse matrices represented in a special format as strings ! :-) Regrettably, all those materials are "missing in action" right now. Those were the days ...

Hope to see you contributing in future S&SM Challenges I may issue, and

Best regards from V.

Hi Raul,

No new mail in the account I told you. Did the message you sent get returned back to you with an "Undelivered" error ?

Please try again. I've tested the address by sending e-mails from other accounts and it does work. If you still experience problems let me know and I'll issue a different address.

Best regards from V.

Hi,

Btw I search actually as much as possible information on the HP8x series (processor, memory controller, ...).

I have an HP85A, an HP86B some 93xx HD and a 9121 floppy. To use the 93xx I need an EMS rom, but can't find one. So I started to dis-assemble and re-assemble roms to allow them to be loaded as binary. I think it's feasible but I need some more information on the 8x series.

For example after searching the web I still have some questions:

- What does SAD and PAD opcodes exactly ? (octal 232 and 237)

- Is there an opcode for octal 326 and 336 ?

- Is there interrupts in the capricorn processor ?

- How the rom and ram paging works in an HP86 ?

- I partially understood the structure of a binary program, but I miss information for the rom structure ...

Thanks in advance if you could help me ,

Olivier

Edited: 1 June 2004, 4:54 a.m.

My server at home must have tipped over; I can't reach it either (for the moment).

Anyway, the listing in directly readable form:

;

;	Assignment #6 for Hewlett-Packard 32SII

;	May 28, 2004

;	was placed in the MoHPC forum

;

LBL Q   ; CK=19C8 032.0

INPUT K

INPUT M

RCL M

1000

STO F

/

RCL M

+

+/-

STO E

FP

ABS

1

+

STO A

LBL A   ; CK=60AE 015.0

RCL K

3

y^x

RCL A

IP

*

STO G   ; A*x^3

RCL E

STO B

LBL B   ; CK=76A0 015.0

RCL K

x^2

RCL B

IP

*       ; B*x^2

RCL+ G

STO H   ; A*x^3+B*x^2

RCL E

STO C

LBL C   ; CK=88E4 013.5

RCL K

RCL C

IP

*       ; C*x

RCL+ H

STO I   ; A*x^3+B*x^2+C*x

RCL E

STO D

LBL D   ; CK=0550 018.0

RCL F

RCL D

IP

RCL+ I  ; A*x^3+B*x^2+C*x+D

x>y?

GTO S   ; remaining D values won't yield

ABS

x<y?

XEQ E   ; more accurate result found

ISG D

GTO D

LBL S   ; CK=5B3E 012.0

ISG C

GTO C

ISG B

GTO B

ISG A

GTO A

STOP

LBL E   ; CK=24F5 018.0

Rolldown

STO F

RCL A

STO P

RCL B

STO Q

RCL C

STO R

RCL D

STO S

RTN

total length:	123.5

;	 storage registers

;	 A B C D parameters

;	 E counting value, typical: -M.00M

;	 F distance to zero of former result

;	 G H I intermediate result to speed up computation

;	 K constant value for x for which the closest root should be found

;	 M upper limit for all of the parameters (-M lower limit for B C D)

;	 P Q R S hold parameters of best result so far

; 

;	 labels

;	 A loop through all A values

;	 B loop through all B values

;	 C loop through all C values

;	 D loop through all D values

;	 E store results of new best result

;	 Q start of program

;	 S to skip irrelevant computations 

Hi, Bram:

Much better, thanks for your considerable efforts to create and share this program, which is a fine example of problem solving using a little HP32-class machine and a lot of user's ingenuity.

Let's hope you'll find my future challenges equally interesting and so you'll be tempted to contribute as well.

Best regards from V.

To recap, the challenge is to:

Quote:

---

"Write a program for your preferred HP calculator to find the integer coefficients A,B,C,D of the cubic equation:

A*x^3 + B*x^2 + C*x + D = 0

which has a root that best approximates a given positive constant K (say, K = Pi), where the coefficients are less or equal in absolute value to a maximum given integer value, M (say, M = 9)."

---

I was looking at the programs offered as solutions for this challenge. While the solutions are certainly clever and very fast, they don't always find the correct polynomial. I don't have all the calculators for which the basic language solutions are provided (HP-71, HP-85, etc.), but I have tested these programs on my Sharp EL-5500III.

The cases originally chosen to test the solutions are:

Quote:

---

Case #1: constant K = Pi, Max. coefficient M = 9, initial interval [3.1, 3.2]
Case #2: same but with Max. coefficient M = 20

---

For K = pi, the solutions presented do seem to return the optimal coefficients. However there are large number of values of K for which the programs do not offer the correct solution.

For example, here are two polynomials for M = 9:

f1(x) = x^3 - 5x^2 + 5x -2, which has a root at x = 3.831177207...

f2(x) = 3x^3 - 9x^2 -8x -6, which has a root at x = 3.832075690...

An exhaustive search shows that there is no other polynomial for the M = 9 case that has a root between these two roots. The midpoint of the roots of these two polynomials is x = 3.831626449...

Now, a correct solution to the challenge would find the coefficients for f1(x) if K is in the range 3.831177207... <= K <= 3.831626449... and the coefficients for f2(x) if K is in the range 3.831626449... < K <= 3.832075690....

However the programs offered as solutions to the challenge find f1(x) for K in the range 3.831177207... <= K <= 3.83192945... and f2(x) for K in the range 3.83192946... <= K <= 3.832075690....

In other words, the incorrect solution is returned for more than 33% of the values of K in the range 3.831177207... <= K <= 3.832075690....

The reason for this is that the programs are trying to find function f(x) that has the minimum value for f(K). This does not necessarily mean that the root for f(x) is closest to K. In this case, the value of K at which f1(K) == f2(K) is not K = 3.831626449... (the midpoint of the two roots), but instead is K = 3.83192945...

I haven't tried the HP-32s solution or the HP-48 solution, but they seem to search for the polynomial in the same way to the Basic language solutions and thus have the same problems. It would therefore appear that the challenge still remains unsolved and that a different approach to an efficient solution is required.

Note that in the original post, Valentin says:

Quote:

---

Your program must input K (say Pi), M (say 9) and a very loose initial interval for the root (say [3.1, 3.2]) and must output the coefficients A,B,C,D of the "best" equation, the value X of the root, and the absolute difference between the root X and the given constant K.

---

Since none of the solutions ask for a loose interval for the root, perhaps that is a clue to obtaining an efficient solution.

Edited: 5 June 2004, 4:00 a.m.

thats most interesting indeed. ive been trying to come up with another way for a submission, but this point throws a spanner there too.

using the same search process, i was going to select vmin, so that whenever f(k) <= vmin, i would run a quick newton solve for f(x)=0 and compare |x-k| minimising this error. the idea was to have vmin small, so that the newton step is quick (eg one iteration). this would give the right answer...

BUT as you've pointed out f(k) could be large and yet |x-k| small with f(x)=0.

so, given its cubic and the M bound, is there is tidy upper bound for vmin that makes this work. naturally, vmin still needs be to quite small otherwise the newton process kicks in on a lot of search cases.

the other question ive been mulling, is whether there's better way altogether than the search process. i had decided to wait for valentin's answer (its less bother if you dont bother).

the only idea i came up with to accelerate the existing search idea is to note that the sign of the coeficients must change at least once when k>0 (eg for pi). this might help eliminate searches when m is larger.

Hi Hugh,

I think you have some good ideas there. One way to get around the problem of large values for f(K), even when f(x) has a root close to K, may be to calculate v = |f(K)/f'(K)| and compare that to the most recent v_min.

There are not a whole lot of extra computations involved, since f'(K) is straightforward to compute (f'(x) = 3ax^2 + 2bx + c). Dividing f(K) by f'(K) should have the effect of 'normalizing' v, making it possible to do the comparison. If f'(x) is a lot less than one, then letting v = f(K) may be fine, avoiding errors that may arise from dividing by a small number.

I'm guessing that it would be necessary to compare |f(K)/f'(K)| to some multiple (2, 3 or more?) of the most recently calculated v_min to ensure that some potential solutions are not eliminated.

It also seems necessary to invoke a root finder each time a suitable polynomial is found to ensure that the correct solution is found. I don't see any easy way around that at the moment. Once a polynomial with a root closer to K has been found, v_min is updated to the v that is calculated for that polynomial.

The fastest solution will eliminate as many of the candidate polynomials as possible before calculating v. Olivier's program does a good job of reducing the number of candidate polynomials, allowing for fast execution times. Your observation of the requirement for a change of sign of at least one of the coefficients could also help execution times.

I probably won't have any time to implement these ideas myself this week, but I'm curious to see what you (or anyone else for that matter) may be able to come up with.

By the way, thanks to Valentin for formulating yet another thought provoking Short & Sweet Math Challenge.

Eamonn.

Hi, Eamonn & Hugh:

Eamon posted:

"By the way, thanks to Valentin for formulating yet another thought provoking Short & Sweet Math Challenge."

Thank you very much, I'm truly glad someone appreciates my humble challenges. Of course it would be great if more people would contribute their ideas, but I guess not so many people as I would desire are keen with things mathematic. Me, it's my #1 lifelong hobby since I can remember.

As for the challenge and the posted solutions you're right, clever and fast as they are they actually minimize the wrong quantity so they'll fail in a large number of instances, as you very cleverly notice. However, the solutions found, while not optimal, have roots very close to the desired value and also very close to the actual optimal solution. In a "production" environment (i.e.: to concoct nice problems for students to solve and such) they would be perfectly ok. But for the theoretical challenge of finding the one optimal solution they don't always deliver.

As stated in my challenge, my original solution was an
HP-71B, 9-line BASIC program (314 bytes), which can solve the two test cases I gave in 7 min. 16 seconds and 3 hour 4 min. respectively. While slower than the posted solutions, mine does actually work in all cases. It goes like this:

  1 DESTROY ALL @ DELAY 0,0 @ DIM M,D,F,X,T,P,C,B,A,K @ INPUT "K,M,P,T=";K,M,P,T

  2 FOR A=1 TO M @ IF SGN(((A*P-M)*P-M)*P-M)=1 THEN END

  3 FOR B=-M TO M @ IF SGN(((A*P+B)*P-M)*P-M)=1 THEN 9

  4 FOR C=-M TO M @ FOR D=-M TO M @ F=SGN(((A*P+B)*P+C)*P+D) @ IF F=SGN(((A*T+B)*T+C)*T+D) THEN 6

  5 X=FNROOT(P,T,((A*FVAR+B)*FVAR+C)*FVAR+D) @ T=ABS(K-X) @ PRINT A;B;C;D;X;T @ P=K-T @ T=K+T @ GOTO 7

  6 IF F=1 THEN 8

  7 NEXT D

  8 NEXT C @ NEXT B

  9 NEXT A

RUN -> K,M,P,T=   

       K    M  P   T               A  B  C  D        X               Err

----------------------------------------------------------------------------

   3.831626,9,3.8,3.9  [INTRO] ->  1 -6  7  5   3.83424318431   .00261718431

                               ->  1 -5  5 -2   3.83117720720   .00044879280

                                           (305 seconds)

   3.831627,9,3.8,3.9  [INTRO] ->  1 -6  7  5   3.83424318431   .00261618431

                               ->  1 -5  5 -2   3.83117720720   .00044979280

                               ->  3 -9 -8 -6   3.83207569042   .00044869042 

                                           (306 seconds)

By the way, Eamonn, congratulations for owning (and obviously liking) such a superb machine as your Sharp EL-5500III. I'm also a Sharp collector and like them a lot. They clearly meet and surpass many HP models past and present, in quality, features and usability as well. May I recommend you also get hold of a Sharp PC-1475, which has all the functionality of your EL-5500III plus a large two-line display and 20-decimal, double precision. Also, the Sharp PC-1350, with its 4-line, fully graphic display is a wonder to behold and use.

Thanks again for your incredibly keen observations and interest and

Best regards from V.

Hi Valentin!

Quote:

---

Of course it would be great if more people would contribute their ideas, but I guess not so many people as I would desire are keen with things mathematic.

---

Actually, I rather suspect that there are many out there like me who appreciate these challenges, and read them, along with the remarkable answers, with great interest---we merely do not have the skills....but we enjoy watching genius unfold, and hope to glean at least a few grains from the experience.

Of course I hope that someday I will manage to have a good attempt that will be worthy of posting!

Best regards,

Bill

Edited: 9 June 2004, 3:30 p.m.

I too, want to add that these challenges are great entertainment. Usually I catch up on the forum while finishing my morning coffee at work. But unfortunately work always presses, and I don't presently have home internet service to pursue the challenge on my own time.

Someday I hope to have a chance to put my solution efforts here, so keep the challenges coming! Thanks, and best regards,

bs

I do really appreciate very much your kind feedback, as it lets me know that my humble efforts do really reach like-minded people, even if they don't choose to post a reply to the particular challenge.

I always aim to be didactic (though at times I may sound 'bombastic' instead) so "readers" are as precious to me as "contributors", to be sure.

More are on their way, stay tuned ! :-)

Best regards from V.

I too, do not always have (or make) time to tackle your challenges, but they are always informative and I do sometimes learn ( admitidly I don't always learn, but thats because many of the solutions are above me).

Nice program to a nifty problem!

Has anybody tried to convert that for the 49 series?

Do you, Valentin, have an example with 10 different values as a result?

[VPN]

Hi Valentin, guys;

I feel I'm a bit of a "math Voyeur", the respectfull kind, of course... >8^O

Please, Valentin, don't take us, at the "silence tribune", as a reference; instead, consider the brave ones that go further and answer your chalenges. When I figure a chance to add some positive comments or possible solutions, I "break the silence" and dare (still in the good sense) posting. Otherwise, I keep my "math Voyerism" in practice...

Best regards from Brazil.

Luiz

Hi, VPN:

VPN posted:

"Do you, Valentin, have an example with 10 different values as a result?"

Excuse me but I don't actually understand your question. Could you elaborate, please ?

P.D.: No disrespect intended but as a Systems Engineer I actually find it funny your VPN moniker (Virtual Private Network] :-)

Never mind!

A) I'm doing the program myself

B) with M=20 one gets 10 "answers"

[Virtual Polish Notation]

I thought last time I saw it your moniker stood for "vertical"(not virtual) polish notation...

identity crisis, anyone?

My true identity revealed:

[Visual Processing Node]

(-;

%%HP: T(3)A(R)F(.);
\<< "SOLVE" { { "K:" "K: guess" 0. } { "P:" "P: min" 0. } { "T:" "T: max" 0. } { "M:" "M: \177loop" 0. } } { } { } { }
IF INFORM
THEN { K P T M } \->TAG DUP LIST\-> 0 DUPDUP TIME \-> K P T M a b I F t
\<< PUSH STD { -1. -55. -72. } SF { -56. -105. } CF CLLCD DATE TIME TSTR 'R.L' DUP2 STO ROT "0" \->TAG 1. \->LIST STO+ 15. 22. SUB DUPDUP 1. 2. SUB { # 0h # 0h } 0. DISPXY ":" { # 7h # 0h } 0. DISPXY 4. 5. SUB { # Ah # 0h } 0. DISPXY ":" { # 11h # 0h } 0. DISPXY 7. 8. SUB { # 14h # 0h } 0. DISPXY "CURR." { # 0h # 7h } 0. DISPXY "-----" { # 0h # Ch } 0. DISPXY "A:" { # 0h # 11h } 0. DISPXY "B:" { # 0h # 18h } 0. DISPXY "C:" { # 0h # 1Fh } 0. DISPXY "D:" { # 0h # 26h } 0. DISPXY 1. M
FOR A " " A DUP 'a' STO + SREV TAIL 1. 3. SUB SREV { # 8h # 11h } 0. DISPXY
IF M P A OVER * PICK3 - OVER * PICK3 - * OVER - SIGN 1. ==
THEN 'A' STO+
ELSE NEG M
FOR B " " B DUP 'b' STO + SREV TAIL 1. 3. SUB SREV { # 8h # 18h } 0. DISPXY
IF M P A OVER * B + OVER * PICK3 - * SWAP - SIGN 1. ==
THEN M 1. + 'B' STO
ELSE M NEG M
FOR C " " C + SREV TAIL 1. 3. SUB SREV { # 8h # 1Fh } 0. DISPXY M NEG M
FOR D " " D + SREV TAIL 1. 3. SUB SREV { # 8h # 26h } 0. DISPXY P A OVER * B + OVER * C + * D + SIGN 'F' STO
IF T A OVER * B + OVER * C + * D + SIGN F ==
THEN
IF F 1. ==
THEN M 1. + 'D' STO
END
ELSE 9999 .1 BEEP K '((A*X+B)*X+C)*X+D' 'X' P T 2. \->LIST ROOT - ABS 'T' STO 'I' INCR 1 - 3 MOD 19 * 9 + R\->B TIME t HMS-
IF DUP 0. <
THEN 24. +
END \-> y h
\<< 'R.L' { A B C D X T h } DUP EVAL 7. \->LIST SWAP \->TAG I \->STR \->TAG 1. \->LIST STO+ "A:" " " A + SREV TAIL 1. 3. SUB SREV + " B:" " " B + SREV TAIL 1. 3. SUB SREV + + " C:" " " C + SREV TAIL 1. 3. SUB SREV + + " D:" " " D + SREV TAIL 1. 3. SUB SREV + + { # 28h } y # 9h - + 0. DISPXY " " I + SREV 1. 2. SUB SREV " X:" + X + { # 1Ch } y # 3h - + 0. DISPXY DATE h TSTR 15. 22. SUB DUPDUP 1. 2. SUB { # 68h } y + 0. DISPXY ":" { # 6Fh } y + 0. DISPXY 4. 5. SUB { # 72h } y + 0. DISPXY ":" { # 79h } y + 0. DISPXY 7. 8. SUB { # 7Ch } y + 0. DISPXY "T:0" T + { # 28h } y # 3h + + 0. DISPXY
\>> K T DUP2 - 'P' STO + 'T' STO
END
NEXT
NEXT
END
NEXT
END
NEXT POP 999. .1 BEEP "=====" { # 0h # 2Ch } 0. DISPXY DATE TIME TSTR 'R.L' OVER STO+ 'R.L' { a b } DUP EVAL 2. \->LIST SWAP \->TAG "\oo" \->TAG 1. \->LIST STO+ 15. 22. SUB DUPDUP 1. 2. SUB { # 0h # 32h } 0. DISPXY ":" { # 7h # 32h } 0. DISPXY 4. 5. SUB { # Ah # 32h } 0. DISPXY ":" { # 11h # 32h } 0. DISPXY 7. 8. SUB { # 14h # 32h } 0. DISPXY 0. WAIT DROP
\>>
END
\>>