 HP50g simplifing a root - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: HP50g simplifing a root (/thread-15660.html) Pages: 1 2 RE: HP50g simplifing a root - Albert Chan - 10-09-2020 05:21 PM I was being stupid. Why solve cubic equation for a This is the way to get cube roots simplified: $$\sqrt{A ± \sqrt{R}} = a ± \sqrt{r} \quad﻿ ⇒ \quad a = \large\frac{\sqrt{A+\sqrt{R}} \;+\; \sqrt{A-\sqrt{R}}}{2}$$ Round a to closest halves, get r, and double check if it round-trip back to (A, R) Code: function simp_cbrt4(A,B,k)  -- simplify cbrt(A + B * sqrt(k))     local R = B*B*k     local a = sqrt(abs(R))     if R<0 then a = (A*A-R)^(1/6) * cos(atan2(a,A)/3)     else        a = (cbrt(A+a) + cbrt(A-a))/2     end     a = a + 0x3p50 - 0x3p50 -- round to halves     local r = (A/a-a*a)/3     local b = B/(3*a*a+r)     if r == b*b*k then return a,b,k end end lua> simp_cbrt4(1859814842094, -59687820010, 415) 11589 ﻿ ﻿ ﻿ ﻿ ﻿ -145 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 415 lua> simp_cbrt4(300940299,103940300,101) 99 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 100 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 101 lua> simp_cbrt(180, 23, 157) ﻿ ﻿ -- (a,b) can be halves 1.5 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.5 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 157 lua> simp_cbrt4(-36, 20, -7) ﻿ ﻿ ﻿ -- work with complex roots too. 3 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 1 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -7 lua> simp_cbrt4(81,30,-3) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -- this is simplest, see comment from previous post. 4.5 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.5 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -3 RE: HP50g simplifing a root - Albert Chan - 10-10-2020 03:58 PM (10-09-2020 05:21 PM)Albert Chan Wrote:  $$\sqrt{A ± \sqrt{R}} = a ± \sqrt{r} \quad﻿ ⇒ \quad a = \large\frac{\sqrt{A+\sqrt{R}} \;+\; \sqrt{A-\sqrt{R}}}{2}$$ Slightly off topics. I was curious how to calculate a accurately (assumed R > 0) Since cube root is odd function, we pull out the sign. a = sign(A)/2 * (³√(|A|+√R) + ³√(|A|-√R)) Last term can be rewritten as ³√(A²-R) / ³√(|A|+√R) Let c = ³√(A²-R), d = ³√(|A|+√R)², we have: a = sign(A)/2 * (d + c) / √d Since both c, d are cube roots, we try identity c³ + d³ = (c + d) (c² - c d + d²) c³ + d³ = (A² - R) + (A² + R + 2|A|√R) = |2A| * (|A|+√R) = |2A| * d√d Substitute back in, many terms get cancelled, a = A / (c²/d - c + d) If ³√(A ± √R) can be simplified, c turns to integer (we had shown this earlier). In fact, the whole denominator turns to integer = a² + 3r = 4a² - 3c Code: function calc_a(A,R)     local c = cbrt(A*A-R)     local d = cbrt(A*A+R + 2*abs(A)*sqrt(R))     return A / (c*c/d - c + d) end lua> function test_a(A,R) : ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ local q = sqrt(R) : ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ return calc_a(A,R), (cbrt(A+q) + cbrt(A-q))/2 : ﻿ ﻿ ﻿ ﻿ ﻿ end lua> lua> test_a(26, 15^2 * 3) 2 ﻿ ﻿ ﻿ ﻿ ﻿ 1.9999999999999982 lua> test_a(9416, 4256^2 * 5) 11 ﻿ ﻿ ﻿ 11.000000000000005 lua> test_a(300940299,103940300^2 * 101) 99 ﻿ ﻿ ﻿ 99.00000000000006 RE: HP50g simplifing a root - Albert Chan - 10-10-2020 04:49 PM (10-10-2020 03:58 PM)Albert Chan Wrote:  Substitute back in, many terms get cancelled, a = A / (c²/d - c + d) If ³√(A ± √R) can be simplified, c turns to integer (we had shown this earlier). In fact, the whole denominator turns to integer = a² + 3r = 4a² - 3c Proof is trivial: if ³√(A ± √R) = a ± √r, then c = ³√(A² - R) = a² - r d = ³√(|A|+√R)² = (|a|+√r)² c²/d + d - c = (|a|-√r)² + (|a|+√r)² - c = (2a² + 2r) - (a² - r) = a² + 3r RE: HP50g simplifing a root - Albert Chan - 10-11-2020 06:28 PM (10-09-2020 02:31 PM)Albert Chan Wrote:  Lets rearrange the cubic to match form x³ + 3px - 2q = 0 c³ = A² - R 4*a³ = 3*c*a + A a³ + 3*(-c/4)*a - 2*(A/8) = 0 → Cubic discriminant = p³ + q² = (-c/4)³ + (A/8)² = (-c³ + A²) / 64 = R/64 → If R < 0, we got 3 real roots. Instead of cubic discriminant, we can show this with more familiar quadratic discriminant. Let f(x) = 4x³ - 3*c*x - A, so that f(a) = 0 f(x) = 4 * (x-a) * (x² + a*x + (a² - 3*c/4)) Quadratic discriminant = a² - 4*(a² - 3*c/4) = -3*(a² - c) = -3*r → if r < 0 (due to R < 0), f got 3 real roots. We can solve the quadratics, but unnecessary. If we have 1 solution, where (a ± √r)³ = A ± √R, we can get the other 2. (a ± √r)³ = ((a ± √r) × w)³ = ((a ± √r) / w)³, ﻿ ﻿ ﻿ ﻿ where w = (-1)^(2/3) Example, simplify ³√(36 + 20i√7) >>> from mpmath import * >>> A, B, k = 36, 20, -7 >>> z = A + B*sqrt(k) >>> q = cbrt(z); print q (3.79128784747792 + 1.27520055582102j) >>> c = cbrt(A*A-B*B*k); print c 16.0 q looks hopeless to simplify, but c is integer, thus possible. Perhaps the "nice" root is not the principle root ? >>> w = exp(2j/3*pi) >>> print q*w, q/w (-3.0 + 2.64575131106459j)﻿ ﻿ ﻿ ﻿ ﻿ ﻿ (-0.791287847477919 - 3.92095186688561j) >>> a = -3 >>> f = lambda x: 4*x**3 - 3*c*x - A >>> print f(a) 0.0 >>> print a, (A/a - a*a)/3 -3 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -7.0 We have ³√(36 + 20i√7) = (-3 + i√7) / (-1)^(2/3) However, converted to form a ± √r, expression look messier than original ³√(A + √R) RE: HP50g simplifing a root - Albert Chan - 10-12-2020 03:17 AM (10-11-2020 06:28 PM)Albert Chan Wrote:  Example, simplify ³√(36 + 20i√7) >>> from mpmath import * >>> A, B, k = 36, 20, -7 ... >>> c = cbrt(A*A-B*B*k); print c 16.0 Here is a novel way to solve for r instead, then get a. (Note: b can be halves, just like a) R/r = (B/b)² = (3c + 4r)² = integer, where r = b²k For this case, we hit the jackpot. If b=±1, 3c+4r = 3c+4k = B r = k = -7 a = A / (c+4r) = 36/ (16-28) = -3 3*arg(a±√r) = 3*atan2(±√7, -3) ≈ ±2.3098 pi We need to match above angle to arg(A+√R), which is on the first quadrant. +2.3098 pi - 2 pi = 0.3098 pi matches. ³√(36 + 20i√7) = (-3 + i√7) / (-1)^(2/3) This is just for fun. No need to solve the cubics. Directly calculate ³√(A+√R) (previous post) is preferred. Comment: Solving for b instead, we can take square root of both side B/b = 3c + 4r = 3c + 4kb² = integer This removed the ambiguity of b = ±1 to b = B/(3c + 4kb²) = 1 Above example need to match angles anyway, so I skipped this. RE: HP50g simplifing a root - peacecalc - 10-12-2020 08:49 PM Hello Albert, I'm really stunnig about your investigations in this topic. Wow, in future I should be more careful in arising such topics! ;-) RE: HP50g simplifing a root - Albert Chan - 10-12-2020 10:54 PM (10-09-2020 02:31 PM)Albert Chan Wrote:  Lets rearrange the cubic to match form x³ + 3px - 2q = 0 c³ = A² - R 4*a³ = 3*c*a + A a³ + 3*(-c/4)*a - 2*(A/8) = 0 → Cubic discriminant = p³ + q² = (-c/4)³ + (A/8)² = (-c³ + A²) / 64 = R/64 → If R < 0, we got 3 real roots. see https://proofwiki.org/wiki/Cardano%27s_Formula I just realized my direct calculation for a is really using Cardano's formula a = ³√(q + √Δ) + ³√(q - √Δ) ﻿ ﻿ ﻿ = ³√((A/8) + √(R/64)) + ³√((A/8) - √(R/64)) ﻿ ﻿ ﻿ = (³√(A+√R) + ³√(A-√R)) / 2 RE: HP50g simplifing a root - CMarangon - 10-12-2020 11:45 PM Hello! Well, I conclude that HP50 does not have the formulas below, into memory. Maybe, it can be included, in next ROM update. Formula: (A ± B⋅√k)^(1/3) = a ± b⋅√k See also post #22, from Albert Chan, above. Carlos (BR) (10-09-2020 02:31 PM)Albert Chan Wrote:   (10-09-2020 12:20 AM)Albert Chan Wrote:  Unfortunately, our divisors based cube root simplify routines were flawed ! Real rational root for a is a divisor of n, *but* possibly divided by 4 ... Another problem is we had assumed there is only 1 real root for a. Lets rearrange the cubic to match form x³ + 3px - 2q = 0 c³ = A² - R 4*a³ = 3*c*a + A a³ + 3*(-c/4)*a - 2*(A/8) = 0 → Cubic discriminant = p³ + q² = (-c/4)³ + (A/8)² = (-c³ + A²) / 64 = R/64 → If R < 0, we got 3 real roots. see https://proofwiki.org/wiki/Cardano%27s_Formula Quote:XCas> find_cbrt(81,30,-3) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → $$\frac{9}{2} + \frac{i}{2}\cdot\sqrt{3}$$ XCas> find_all_a(A,R) := solve(surd(A*A-R,3) = a*a - (A/a-a*a)/3 , a) XCas> find_all_a(81, 30^2 * -3) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → [-3, -3/2, 9/2] XCas> simplify( [-3+2i*sqrt(3), -3/2-5i/2*sqrt(3), 9/2+i/2*sqrt(3)] .^ 3) $$\qquad\quad [81 + 30i\sqrt3\;,\; 81 + 30i\sqrt3\;,\; 81 + 30i\sqrt3$$ If we consider integer as simplest, then $$\sqrt{81 + 30i\sqrt3} = -3+2i\sqrt3$$ Comment: I was wrong on above example. x³ = y does not imply ³√y = x. We need to consider branch-cut In other words, simplified form should match numerical evaluated values. >>> (81 + 30j * 3**0.5) ** (1/3) (4.5+0.86602540378443837j) RE: HP50g simplifing a root - peacecalc - 10-13-2020 06:30 AM Hello Albert, Quote:a = ³√(q + √Δ) + ³√(q - √Δ) ﻿ ﻿ ﻿ = ³√((A/8) + √(R/64)) + ³√((A/8) - √(R/64)) ﻿ ﻿ ﻿ = (³√(A+√R) + ³√(A-√R)) / 2 This was the origin of my thoughts for simplyfing, because I construct examples to test cardanos formular i keyed in my HP50g and I got that kind of "garbage" for simple solutions i. e. 1 minus root of two. So my question arises, why my calc cannot simplify this...the rest is history in nearly 30 posts... RE: HP50g simplifing a root - peacecalc - 10-13-2020 06:36 AM Hello CMaragon, don't wait for a rom update, I'm feared you don't live long enough, because the HP50g is no longer produced and so the last rom version is 2.15. It's a pity, but for HP calc department is "prime time". RE: HP50g simplifing a root - grsbanks - 10-13-2020 06:46 AM (10-12-2020 11:45 PM)CMarangon Wrote:  Well, I conclude that HP50 does not have the formulas below, into memory. Maybe, it can be included, in next ROM update. Sadly the 50g was discontinued 5 years ago. There will be no ROM update RE: HP50g simplifing a root - Albert Chan - 10-24-2020 02:19 PM (10-11-2020 06:28 PM)Albert Chan Wrote:  We have ³√(36 + 20i√7) = (-3 + i√7) / (-1)^(2/3) Mathematica gives many "simplifed" forms, here are 2 of them. (-3+i√7) * (-1-i√3)/2 = (√21+3)/2 + i/2*(3√3-√7) = (√21+3)/2 + i/2*√(34-6√21) I was curious, given √(34-6√21), how to "simplify" to (3√3-√7) ? (10-09-2020 05:21 PM)Albert Chan Wrote:  $$\sqrt{A ± \sqrt{R}} = a ± \sqrt{r} \quad﻿ ⇒ \quad a = \large\frac{\sqrt{A+\sqrt{R}} \;+\; \sqrt{A-\sqrt{R}}}{2}$$ We can use the same trick (if LHS is real, for both ±) $$\sqrt{A ± B\sqrt{c d}} = a\sqrt{c} \,±\, b\sqrt{d} \quad﻿ ⇒ \quad a\sqrt{c} = \large\frac{\sqrt{A+B\sqrt{c d}}\;+\;\sqrt{A-B\sqrt{c d}}}{2}$$ lua> A, B, k = 34, -6, 21 lua> C = B * sqrt(k) lua> mean = (sqrt(A+C) + sqrt(A-C))/2 lua> mean^2 27 27 = 3*3*3 = a*a*c. We try a=c=3, d=7, so just solve for b: (3√3 + b√7)² = (27+7b²) + 6b√21 = 34 - 6√21 We have b = -1, thus √(34 - 6√21) = 3√3 - √7 It is nice that XCas do this automatically XCas> simplify(√(34-6*√(21))) ﻿ ﻿ ﻿ ﻿ ﻿ → 3√3 - √7