Quadratic Solver for Casio FX-3650P

[Albert Chan]

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

[Albert Chan]

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

[Albert Chan]

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

[Albert Chan]

(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

[Albert Chan]

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

[Albert Chan]

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

[Albert Chan]

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

[Albert Chan]

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

[Albert Chan]

(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

[Albert Chan]

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