Quadratic Solver for Casio FX-3650P - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: Not HP Calculators (/forum-7.html) +--- Forum: Not remotely HP Calculators (/forum-9.html) +--- Thread: Quadratic Solver for Casio FX-3650P (/thread-11075.html) Quadratic Solver for Casio FX-3650P - Albert Chan - 07-17-2018 04:29 PM 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 ... RE: Quadratic Solver for Casio FX-3650P - Albert Chan - 07-17-2018 10:52 PM 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``` RE: Quadratic Solver for Casio FX-3650P - Albert Chan - 07-18-2018 04:35 PM 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 RE: Quadratic Solver for Casio FX-3650P - Albert Chan - 08-10-2018 04:28 PM (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 RE: Quadratic Solver for Casio FX-3650P - Albert Chan - 08-11-2018 02:29 PM 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 RE: Quadratic Solver for Casio FX-3650P - Albert Chan - 08-11-2018 10:13 PM 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 RE: Quadratic Solver for Casio FX-3650P - Albert Chan - 11-09-2018 06:43 PM 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 RE: Quadratic Solver for Casio FX-3650P - Albert Chan - 11-10-2018 08:14 PM 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 RE: Quadratic Solver for Casio FX-3650P - Albert Chan - 11-11-2018 05:44 PM (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 RE: Quadratic Solver for Casio FX-3650P - Albert Chan - 12-03-2018 12:08 PM 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