Threaded Mode | Linear Mode

Roberto Volpi · 05-21-2023, 10:49 AM

(05-20-2023 12:05 PM)Thomas Klemm Wrote: It took me a while to understand how your program works.
The following Python program helped me to do so:

Code:

from math import copysign def cbeq(a, b, c, d): B, C, D = b / a, c / a, d / a # C002 - C011 F = 0 H = copysign(1, D) # C012 - C014 X = abs(D) + abs(B) # C015 - C019 P = 10 # C020 - C022 O = 1e-6 # C023 - C024 Y, X = X, O X *= P # C025 print(f"C025: {X=:15.12f} {Y=:15.12f} {X <= Y}") while X <= Y: # C026 - C027 X *= P # C025 print(f"C025: {X=:15.12f} {Y=:15.12f} {X <= Y}") K = X # C028 L = 0 while True: if L * H >= 0: # C054 - C058 K /= P # C029 - C030 H = -H # C031 - C033 G = K * H + F # C034 - C037 if G == F: # C038 - C039 return G # C040 L = ((G + B) * G + C) * G + D # C041 - C047 print(f"C048: {H=:2.0f} {K=:15.12f} {G=:15.12f} {L=:15.12f}") if O > abs(L): # C048 - C050 return G # C051 F = G # C052 - C053
First it tries to find an integral power of 10 that is bigger or equal to \(|\frac{B}{A}| + |\frac{D}{A}|\).
Then a digit by digit algorithm is used to approach the solution.
Whenever a change of sign is found the next decimal place to the right is considered.

Here's an example for the equation \(x^3 - x^2 - x - 1 = 0\):

cbeq(1, -1, -1, -1)

Code:

C025: X= 0.000010000000 Y= 2.000000000000 True C025: X= 0.000100000000 Y= 2.000000000000 True C025: X= 0.001000000000 Y= 2.000000000000 True C025: X= 0.010000000000 Y= 2.000000000000 True C025: X= 0.100000000000 Y= 2.000000000000 True C025: X= 1.000000000000 Y= 2.000000000000 True C025: X=10.000000000000 Y= 2.000000000000 False C048: H= 1 K= 1.000000000000 G= 1.000000000000 L=-2.000000000000 C048: H= 1 K= 1.000000000000 G= 2.000000000000 L= 1.000000000000 C048: H=-1 K= 0.100000000000 G= 1.900000000000 L= 0.349000000000 C048: H=-1 K= 0.100000000000 G= 1.800000000000 L=-0.208000000000 C048: H= 1 K= 0.010000000000 G= 1.810000000000 L=-0.156359000000 C048: H= 1 K= 0.010000000000 G= 1.820000000000 L=-0.103832000000 C048: H= 1 K= 0.010000000000 G= 1.830000000000 L=-0.050413000000 C048: H= 1 K= 0.010000000000 G= 1.840000000000 L= 0.003904000000 C048: H=-1 K= 0.001000000000 G= 1.839000000000 L=-0.001568281000 C048: H= 1 K= 0.000100000000 G= 1.839100000000 L=-0.001021459529 C048: H= 1 K= 0.000100000000 G= 1.839200000000 L=-0.000474547712 C048: H= 1 K= 0.000100000000 G= 1.839300000000 L= 0.000072454457 C048: H=-1 K= 0.000010000000 G= 1.839290000000 L= 0.000017750174 C048: H=-1 K= 0.000010000000 G= 1.839280000000 L=-0.000036953205 C048: H= 1 K= 0.000001000000 G= 1.839281000000 L=-0.000031482908 C048: H= 1 K= 0.000001000000 G= 1.839282000000 L=-0.000026012602 C048: H= 1 K= 0.000001000000 G= 1.839283000000 L=-0.000020542286 C048: H= 1 K= 0.000001000000 G= 1.839284000000 L=-0.000015071962 C048: H= 1 K= 0.000001000000 G= 1.839285000000 L=-0.000009601629 C048: H= 1 K= 0.000001000000 G= 1.839286000000 L=-0.000004131286 C048: H= 1 K= 0.000001000000 G= 1.839287000000 L= 0.000001339065 C048: H=-1 K= 0.000000100000 G= 1.839286900000 L= 0.000000792030
And we end up with:

G= 1.839286900000

But this made me think, what happens if the solution gets closer and closer to 1e-6.

cbeq(1e3, -1e2, -1e1, -1)

Code:

C025: X= 0.000010000000 Y= 0.101000000000 True C025: X= 0.000100000000 Y= 0.101000000000 True C025: X= 0.001000000000 Y= 0.101000000000 True C025: X= 0.010000000000 Y= 0.101000000000 True C025: X= 0.100000000000 Y= 0.101000000000 True C025: X= 1.000000000000 Y= 0.101000000000 False C048: H= 1 K= 0.100000000000 G= 0.100000000000 L=-0.002000000000 C048: H= 1 K= 0.100000000000 G= 0.200000000000 L= 0.001000000000 C048: H=-1 K= 0.010000000000 G= 0.190000000000 L= 0.000349000000 C048: H=-1 K= 0.010000000000 G= 0.180000000000 L=-0.000208000000 C048: H= 1 K= 0.001000000000 G= 0.181000000000 L=-0.000156359000 C048: H= 1 K= 0.001000000000 G= 0.182000000000 L=-0.000103832000 C048: H= 1 K= 0.001000000000 G= 0.183000000000 L=-0.000050413000 C048: H= 1 K= 0.001000000000 G= 0.184000000000 L= 0.000003904000 C048: H=-1 K= 0.000100000000 G= 0.183900000000 L=-0.000001568281 C048: H= 1 K= 0.000010000000 G= 0.183910000000 L=-0.000001021460 C048: H= 1 K= 0.000010000000 G= 0.183920000000 L=-0.000000474548

G= 0.183920000000

cbeq(1e6, -1e4, -1e2, -1)

Code:

C025: X= 0.000010000000 Y= 0.010001000000 True C025: X= 0.000100000000 Y= 0.010001000000 True C025: X= 0.001000000000 Y= 0.010001000000 True C025: X= 0.010000000000 Y= 0.010001000000 True C025: X= 0.100000000000 Y= 0.010001000000 False C048: H= 1 K= 0.010000000000 G= 0.010000000000 L=-0.000002000000 C048: H= 1 K= 0.010000000000 G= 0.020000000000 L= 0.000001000000 C048: H=-1 K= 0.001000000000 G= 0.019000000000 L= 0.000000349000

G= 0.019000000000

cbeq(1e9, -1e6, -1e3, -1)

Code:

C025: X= 0.000010000000 Y= 0.001000001000 True C025: X= 0.000100000000 Y= 0.001000001000 True C025: X= 0.001000000000 Y= 0.001000001000 True C025: X= 0.010000000000 Y= 0.001000001000 False C048: H= 1 K= 0.001000000000 G= 0.001000000000 L=-0.000000002000

G= 0.001000000000

Unfortunately it looks like fewer and fewer digits are calculated.

Long story short, the algorithm to calculate the first, always real root, is an open numerical method.

We start with a guess, and decimals tlll 10E-6 are calculated.

I have also tried to use native solver of HP35S. The program was shorter, but sometimes the accuracy was poor, sometime execution time was too long, or both.

I have used the constant value 1E-6 in order to balance speed and accuracy. I consider extremely important to obtain the first root in a decent time, but those who prefer a different value may change instruction C023 as they like.