HP Forums - Quadratic Solver for Casio FX-3650P

Quadratic Solver had been done before, but I wanted accurate roots.

The school taught (-b +/- SQRT(b^2 - 4ac)) / (2a) may have one root that loses precision.

Example: z^2 - (5e9 + 1) z + 5e9 should have roots = 5e9 and 1
But, on my FX-115MS solver, I get roots = 5e9, 1.001 (only 3 digits accurate)

Besides, my FX-3650P does not have a quadratic solver, so here it is: (45 keys)

-- Solve for z: A z^2 + B z + C = 0
-- If real roots, roots saved in X, Y
-- If MATH ERROR (SQRT of negative number), complex roots = X(1 +/- SQRT(Y))

Code:

? -> A : ? -> B : ? -> C :

-B / 2A -> X :

1 - 4AC / B^2 -> Y :

X + X SQRT Y -> X PAUSE

C / AX -> Y

Redo above example, I get roots = 5e9, 1 :-)

Besides accuracy, other benefits of this program:

1. Roots are saved in variable X, Y, keeping internal precision.
-- A, B, C are also saved

2. Because of better precision, complex roots easier to recover its radical form:
-- Example roots of 3 z^2 + 7 z + 10 = 0:

Roots (from above program, recover fraction with FRAC key)
= (-1.166666667)(1 +/- SQRT(-1.448979592))
= (-7/6)(1 +/- SQRT(-71/49))
= -7/6 +/- SQRT(71)/6 i

3. Complex roots is warned with MATH ERROR (sqrt of negative number), roots = X(1 +/- SQRT(Y))

FX-115MS Solver only showed a tiny dot (r<=>i) on the upper right screen.
I only noticed it when both roots look exactly the "same"
Finally, I realized the roots are complex, but have the same real part ...

Looks like solving quadratic equation accurately is harder than it look ...

Thomas Klemm article about solving quadratic equation is educational
http://www.hpmuseum.org/forum/thread-2181.html

Because FX-3650P does not have a ABS function, I had rearranged equation to produce the "bigger" root.

(-B +/- SQRT(B^2 - 4AC)) / (2A) = -B/(2A) * (1 +/- SQRT(1 - 4AC/B^2))

Bigger root (ignore sign) = -B/(2A) (1 + SQRT(1 - 4AC/B^2)) = X(1 + SQRT Y)

The problem with above is it may have worse precision for the "discriminant" Y

This is my second attempt (86 keys)

Code:

? -> A : ? -> B : ? -> C : B^2 - 4AC -> D : 

D >= 0 => GOTO 1 : -- real roots branch

-B/2A -> X PAUSE -- this version, return X +/- Y i for complex roots

SQRT -D / 2A -> Y PAUSE  -- Fall thru LBL1 to quit with MATH ERROR

LBL 1: SQRT D -> Y

B >= 0 => GOTO 2 : -Y -> Y : -- make sure B and Y have same sign

LBL2: -(B + Y) / 2A -> X PAUSE -- big root

C / AX -> Y -- small root

Using Thomas Klemm idea about coefficient reduction, my Quadratic Solver can now handle harder Quadratic.
http://www.hpmuseum.org/forum/thread-2181.html

For quadratic with a hard time getting good value of discriminant B^2 - 4AC
Transform Quadratic F(z) with offset k, so z = y + k:

F(z) = A z^2 + B z + C
= A(y^2 + 2ky + k^2) + B(y + k) + C
= A y^2 + (2kA + B) y + (A k^2 + B k + C)
= A y^2 + F'(k) y + F(k)
= G(y)

To reduce the G(y) coefficient F'(k) ~ 0, k ~ -B/(2A)
After solving G(y) = 0 for y, z = y + k, and we are done

To automate above, I made another program for F(z) / (z - X)
Each coefficient entry showed the coefficient of F(z) / (z - X)
After finished entering the coefficients, B = F'(X), C = F(X)

Since the program calculated F'(X) and F(X), it also does Newtons method.
If X=0, the program assumed I wanted Newton improved X, based on previous results.

Synthetic Division program listing (59 keys):

Code:

? -> X : X=0 => Y - C/B -> X PAUSE -- Newton correction

X -> Y : 0 -> B : ? -> A : A -> D :

LBL 1 : D -> C : ? -> C : -- Horner's Method: F(z) / (z - X)

BX + D -> B : -- F'(X)

DX + C -> D : -- F(X)

GOTO 1 :       -- Key AC to stop

Borrowing Klemm example, solve F(z) = 11713 z^2 - 1470492 z + 46152709 = 0

Just by looking at the coefficients, my Casio cannot give exact B^2 - 4AC
So, try reduction with offset X ~ -B/(2A) = 62.77179203

We wanted the *exact* transformed G(y), so pick X = 63
Running above program, we get:

Synthetic Division Program:
? X -- 63 enter
? A -- 11713 enter
? C -- -1470492 enter
? C -- 46152709 enter AC (to stop)

--> A = 11713, B = 5346, C = 610
--> Now, solve G(y) = 11713 y^2 + 5346 y + 610 = 0

Quadratic Solver Program:
? A -- enter
? B -- enter
? C -- enter
MATH ERROR -- complex roots = X +/- Y i

--> G(y) roots = -0.228207974 +/- 8.537522411e-05 i
--> F(z) roots = 62.771792026 +/- 8.537522411e-05 i

(07-18-2018 04:35 PM)Albert Chan Wrote: [ -> ]Borrowing Klemm example, solve F(z) = 11713 z^2 - 1470492 z + 46152709 = 0

Just by looking at the coefficients, my Casio cannot give exact B^2 - 4AC
So, try reduction with offset X ~ -B/(2A) = 62.77179203 ...

With a little bit of work, above discriminant can be calculated exactly
Say, a calculator only have 10 digits precision, so any 5-digits multiply are exact:

B^2 = 1470492^2
= (14e5 + 70492)^2
= (14e5 + 2*70492) * 14e5 + 4969122064
= (21573776 + 49691) * 1e5 + 22064
= 21623467e5 + 22064

4 A C = 4 * 11713 * 46152709
= 46852 * (461e5 + 52709)
= 21598772e5 + 2469522068
= 21623467e5 + 22068

D = B^2 - 4 A C = 22064 - 22068 = -4

I shown the steps, but on the calculator, it is much simpler:
B^2 - 4AC = (14e5 + 70492)^2 - 46852 * (461e5 + 52709)

1e5 terms:
(14e5 + 2*70492)*14e5 - 46852*461e5 => -24996e5

D = (reduced B^2 - 4AC) + (1e5 terms):
70492^2 - 46852*52709 + Ans => -4

Here is a revised Quadratic Solver (85 steps), discriminant D adjustable (if better available)
Since D is shown, it's sign signal real or complex roots.

Code:

? -> A : ? -> B : ? -> C : B^2 - 4 A C -> D : 

? -> D : D >= 0 => GOTO 1 : 

-B / 2A -> X PAUSE SQRT -D / 2A -> Y PAUSE

LBL 1 : SQRT D -> Y : B > 0 => -Y -> Y : Y - B -> Y : 

2C / Y -> X PAUSE

Y / 2A -> Y

With this, and bits of work to get exact D, I can solve Cadillac Quadratic Solver, case 8
(Casio FX-3650P don't have enough precision to handle, even with reduced B, C)

Get D with previous post 1e5 trick
B^2 - 4 A C = (-22222222)^2 - 4 * 8441600 * 14624809
= (222e5 + 22222)^2 - (337e5 + 66400) (146e5 + 24809)

1e5 terms:
(222e5 + 2*22222) * 222e5 - 337e5 * 146e5 - 337e5 * 24809 - 66400 * 146e5 => 11535e5

D = (reduced B^2 - 4 A C) + (1e5 terms)
(22222^2 - 66400 * 24809) + ANS => -316

Run above:
? A -- Enter 8441600
? B -- Enter -22222222
? C -- Enter 14624809
? D -- Enter -316 (calculated D were bad, shown -1000)

-> Roots = 1.316232823 +/- 1.052904001e-6 * I

There is a faster way to figure out exact value of B^2 and 4AC, thus D
(No need to collect all those 1e5 terms)

D = B^2 - 4 A C
= (-22222222)^2 - 4 * 8441600 * 14624809
= 22222222^2 - 33766400 * 14624809

On the Casio, both B^2 and 4AC get the same result R = 4938271506e5

B^2 % 1e5 = 22222^2 % 1e5 = 17284
4AC % 1e5 = 66400*24809 % 1e5 = 17600

D = B^2 - 4AC = (R + 17284) - (R + 17600) = -316

Another way to get discriminant, from Kahans cubic equation paper (page 11)
BTW, page 12 example 2 typo, should be: 1,160,928,203² - 3,234,424,451 * 416,690,636 = -89,060,331,627

PHP Code:

function b2_less_ac(a, b, c)  -- assume a, c > 0
    if a < c then a, c = c, a end
    if verbal then io.write(math.abs(b), '\xfd - ', a, ' * ', c, '\n') end
    local n = 0
    if c>0 then n = math.floor(b/c + 0.5) end
    if n == 0 then return b*b - a*c end
    local h = math.floor(n/2)
    a = a - (n-h)*b - h*b
    b = b - (n-h)*c - h*c
    return b2_less_ac(c, b, a - n*b)
end 

To show that reduction code work, let reduced a = a', reduced b = b':
b'² = (b - nc)² = b² - nc * (2b - nc) = b² - nc * (b+b')
a'c = (a - nb - nb') * c = ac - nc * (b+b')
-> b'² - a'c = b² - ac

lua> A, B, C = 8441600, -22222222, 14624809 -- previous post example
lua> verbal = true
lua> b2_less_ac(4*A, B, C)
22222222² - 33766400 * 14624809
7027396² - 14624809 * 3376748
273900² - 3376748 * 22217
7296² - 22217 * 2396
108² - 2396 * 5
2² - 64 * 5
-316

Doing by hand, we can do slightly better:

(-22222222)² - (4)(8441600)(14624809)
= 22222222² - (5849923600)(84416) -- n = round(22222222/84416) = 263, b' = b-nc = 20814
= 20814² - (84416)(5849923600 - 263*22222222 - 263*20814)
= 20814² - (84416)(5132)
= -316

Discovered a trivia: discriminant is the same if quadratic is "shifted"

AX² + BX + C, let Y = X - k

Using synthetic division, we get: AY² + B'Y + C', where B' = 2Ak + B, C' = Ak² + Bk + C

B'² - 4 A C'
= (2Ak + B)² - 4A*(Ak² + Bk + C)
= (4A²k² + 4ABk + B²) - (4A²k² + 4ABk + 4AC)
= B² - 4 A C

Using previous post example: A, B, C = 8441600, -22222222, 14624809

-B/(2A) ≈ 1.316232823
Let Y = X - 1.3, quadratic => 8441600Y² - 274062 Y + 2224.4

4AC' has 10 significant digits, thus can evaluate exactly.
Slight adjustment made B'² part also exact ...

B'² - 4 A C'
= -4 A C' + (B'-2)(B'+2) + 4
= (-4)(8441600)(2224.4) + (274060)(274064) + 4
= -316

Edit: there is a simpler prove that discriminant unchanged when shifted.
Shifting is just another perspective to say where is considered zero.
The gap between the roots remains the same.
Since gap = √(D / A²), and A unchanged when shifted, D also unchanged.

Update: we could avoid B' adjustment by doing synthetic division again
let Z = Y - 0.016, quadratic => 8441600Z² - 3930.8 Z + 0.4576 = 0
Solve for Z, then X = Z + 1.316

(11-10-2018 08:14 PM)Albert Chan Wrote: [ -> ]Discovered a trivia: discriminant is the same if quadratic is "shifted"

AX² + BX + C, let Y = X - k

Using synthetic division, we get: AY² + B'Y + C', where B' = 2Ak + B, C' = Ak² + Bk + C

Using this trivia, calculator can get a more precise √2, positive root of X² - 2
Start with k ~ √2 = 1.414213562, Y = X - k

A = 1
B' = 2k = 2.828427124
C' = k² - 2 ~ -1.06e-9, but we need this exact ...
D = B² - 4 A C = 0 - 4(1)(-2) = 8

C' = (1.4142 + 13562e-9)² - 2
= -2 + 1.4142² + 2*1.4142*13562e-9 + (13562e-9)²
= -1.055272156e-9

Another approach is reduce B² and AC, a bit at a time:

1e18 C'
= 1414213562² - (1e10)(2e8) -- n = 7, b-nc = 14213562
= 14213562² - (2e8)(1e10 - 7*(1414213562+14213562))
= 14213562² - (2e8)(1010132)
= 14213562² - (2e3)(1010132)(1e5) -- n = 142, b-nc = 13562
= 13562² - (1e5)(2020264e3 - 142*(14213562 + 13562))
= -1055272156

Modulo math, using Chinese Remainder Theorem, also work well:

Let a ≡ 1e18 C' (mod 1e5) ≡ 13562² ≡ 27844
Let b ≡ 1e18 C' (mod 1e5 + 1) ≡ (-14142+13562)² - (1)(-2e3) ≡ 38397
→ 1e18 C' (mod 1e10 + 1e5) ≡ (a-b)*1e5 + a ≡ -1055272156

Enter A, B', C', and D to the Solver, we get Y correction = 3.730950488e-10

√2 ~ 1.414213562 3730950488

Just learn about correction for CCDF function, perhaps useful for above √(2) correction:

Redo above, without solving the quadratic: 0 = y² + 2.828427124 y - 1.055272156e-9 = a y² + b y + c

1st order correction: y ~ t = -c/b
Casio (with internal precision digits), y = 3.73095048851e-10

2nd order correction: y ~ t - (a/b) t²
Casio (with internal precision digits), y = 3.73095048802e-10

√2 ~ 1.414213562 3730950488 02

Edit: Assume a y² + b y + c small real root exist, t = -c/b, w = -(a/b) t = ac/b² :

y = t * (1 + w + 2w^2 + 5w^3 + 14w^4 + 42w^5 + 132w^6 + ...)

Redo above with 40 digits precision: y1 = t, y2 = t + tw ...

y1: 3.730950488509033262969090378444553482510E-10
y2: 3.730950488016887241967143787508642382276E-10
y3: 3.730950488016887242096980785739535431605E-10
y4: 3.730950488016887242096980785696718753754E-10