[PATCH] 2/2 simplify square root in 'radicalSolve'

oldk1331

unread,

Apr 8, 2019, 8:23:54 AM4/8/19

to fricas...@googlegroups.com

This patch simplifies the root of "x^2+2*a*x+2*b=0":

(1) -> rhs first radicalSolve(x^2+2*a*x+2*b,x)

+------------+
| 2
- \|- 8 b + 4 a - 2 a
(1) -----------------------
2
Type: Expression(Integer)
to

(2) -> -sqrt(a^2-2*b)-a

+----------+
| 2
(2) - \|- 2 b + a - a
Type: Expression(Integer)

by using sqaure-free factorization of Expression.

https://github.com/oldk1331/fricas/commit/650e788da9c5247658db6115f80bc068f68ff870.patch

diff --git a/ChangeLog b/ChangeLog
index 160f1cf9..66bea43d 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2019-04-08 Qian Yun <oldk...@gmail.com>
+
+ * src/algebra/solvefor.spad, src/algebra/solverad.spad:
+ simplify square root in 'radicalSolve'
+
2019-04-08 Qian Yun <oldk...@gmail.com>

* src/algebra/solvefor.spad: remove unused 'mapSolve'
diff --git a/src/algebra/solvefor.spad b/src/algebra/solvefor.spad
index 2bab4943..dab31527 100644
--- a/src/algebra/solvefor.spad
+++ b/src/algebra/solvefor.spad
@@ -1,6 +1,8 @@
)abbrev package SOLVEFOR PolynomialSolveByFormulas
--- Examples of fields with "^": (%, Fraction Integer) -> % are
--- Complex Float, RealClosure(K) and AlgebraicNumber
+-- Examples of \spadtype{F} (fields with "^" : (%, Fraction Integer) -> %) are
+-- Expression Integer, Complex Float, RealClosure(K) and AlgebraicNumber
+-- If \spadtype{F} is Expression(RR), you should put \spadtype{RR} as
+-- the third argument to this package for simplification purposes.
-- RealClosure(K) is unlikly to work here...
++ Author: Stephen M. Watt, Barry Trager
++ Description:
@@ -9,10 +11,11 @@
++ Care is taken to introduce few radicals so that radical extension
++ domains can more easily simplify the results.

-PolynomialSolveByFormulas(UP, F) : PSFcat == PSFdef where
+PolynomialSolveByFormulas(UP, F, RR) : PSFcat == PSFdef where

UP : UnivariatePolynomialCategory F
F : Field with "^": (%, Fraction Integer) -> %
+ RR : Join(PolynomialFactorizationExplicit, Comparable, CharacteristicZero)

L ==> List

@@ -55,6 +58,23 @@
part(s : F) : F ==
s

+ -- if F is Expression(RR), then we can simplify "s^(1/2)" to
+ -- "a*b^(1/2)" by sqaure-free factorization. Note that we should
+ -- the unit separately; and the sign of "a" is not important --
+ -- we will return both roots later.
+ sqrt2(s : F) : F ==
+ if F is Expression RR then
+ smpsqfr := squareFree numer s
+ a : F := coerce("*"/[f.factor^(f.exponent quo 2) for f in factorList smpsqfr])
+ b : F := coerce("*"/[f.factor^(f.exponent rem 2) for f in factorList smpsqfr])
+ u : RR := retract(coerce(unit smpsqfr)@F)
+ sqfru := squareFree u
+ ua := "*"/[f.factor^(f.exponent quo 2) for f in factorList sqfru]
+ ub := "*"/[f.factor^(f.exponent rem 2) for f in factorList sqfru]
+ ua*a*sqrt(ub*b)
+ else
+ s^(1/2)
+
-----------------------------------------------------------------
-- Entry points and error handling
-----------------------------------------------------------------
@@ -139,7 +159,7 @@
needLcoef c2; needChar0()
(c0 = 0) => cons(0$F, linear(c2, c1))
(c1 = 0) => [(-c0/c2)^(1/2), -(-c0/c2)^(1/2)]
- D := part(c1^2 - 4*c2*c0)^(1/2)
+ D := sqrt2(c1^2 - 4*c2*c0)
[(-c1+D)/(2*c2), (-c1-D)/(2*c2)]

aQuadratic(c2, c1, c0) ==
diff --git a/src/algebra/solverad.spad b/src/algebra/solverad.spad
index 46a5cf7b..50fba034 100644
--- a/src/algebra/solverad.spad
+++ b/src/algebra/solverad.spad
@@ -76,7 +76,7 @@
L ==> List
P ==> Polynomial

- SOLVEFOR ==> PolynomialSolveByFormulas(SUP RE, RE)
+ SOLVEFOR ==> PolynomialSolveByFormulas(SUP RE, RE, R)
UPF2 ==> SparseUnivariatePolynomialFunctions2(PR, RE)

Cat ==> with
diff --git a/src/input/bugs2019.input b/src/input/bugs2019.input
index f3268439..db664c40 100644
--- a/src/input/bugs2019.input
+++ b/src/input/bugs2019.input
@@ -20,4 +20,9 @@ testEquals("factor(x-1)*0", "0")
testEquals("gcd(factor(x-1), 1)", "1")
testEquals("factor(x^2-1)+1-x^2", "0")

+testcase "simplification of square root in 'radicalSolve'"
+
+testEquals("rhs first radicalSolve(x^2+2*a*x+2*b,x)", "-sqrt(a^2-2*b)-a")
+testEquals("rhs first radicalSolve(x^2+2*a*c*x+2*b*c^2,x)", "-c*sqrt(a^2-2*b)-a*c")
+
statistics()

oldk1331

unread,

Apr 13, 2019, 9:06:30 PM4/13/19

to fricas...@googlegroups.com

In "sqrt2", instead of factoring expression manually, we can use existing "froot":

--- a/src/algebra/solvefor.spad
+++ b/src/algebra/solvefor.spad
@@ -64,14 +64,10 @@

-- we will return both roots later.

sqrt2(s : F) : F ==

if F is Expression RR then

- smpsqfr := squareFree numer s
- a : F := coerce("*"/[f.factor^(f.exponent quo 2) for f in factorList smpsqfr])
- b : F := coerce("*"/[f.factor^(f.exponent rem 2) for f in factorList smpsqfr])
- u : RR := retract(coerce(unit smpsqfr)@F)
- sqfru := squareFree u
- ua := "*"/[f.factor^(f.exponent quo 2) for f in factorList sqfru]
- ub := "*"/[f.factor^(f.exponent rem 2) for f in factorList sqfru]
- ua*a*sqrt(ub*b)
+ K ==> Kernel F
+ PP ==> SparseMultivariatePolynomial(RR, K)
+ res := froot(s, 2)$PolynomialRoots(IndexedExponents K, K, RR, PP, F)
+ res.coef * (res.radicand)^(1/res.exponent)
else
s^(1/2)

Waldek Hebisch

unread,

Apr 16, 2019, 2:41:44 PM4/16/19

to fricas...@googlegroups.com

oldk1331 wrote:
>
> In "sqrt2", instead of factoring expression manually, we can use existing
> "froot":

OK. But this and previous patch really should go as one. So please
commit as a single change.

> --
> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.
> To post to this group, send email to fricas...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/CAGBJN93gf%3Du2J8MkDi3QPvNC7wVn2BHAfEnXUcM6-UHM_7DZWA%40mail.gmail.com.
> For more options, visit https://groups.google.com/d/optout.
>
> --0000000000004e2b0f05867325e0
> Content-Type: text/html; charset="UTF-8"
> Content-Transfer-Encoding: quoted-printable
>
> <div dir=3D"auto">In "sqrt2", instead of factoring expression man=
> ually, we can use existing "froot": 
> 
> --- a/src/algebra/solvefor.spad 
> +++ b/src/algebra/solvefor.spad 
> @@ -64,14 +64,10 @@ 
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0-- we will return both roots later. 
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0sqrt2(s : F) : F =3D=3D 
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if F is Expression RR then<=
> br>
> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 smpsqfr :=3D squar=
> eFree numer s 
> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 a : F :=3D coerce(=
> "*"/[f.factor^(f.exponent quo 2) for f in factorList smpsqfr])<br=
> >
> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 b : F :=3D coerce(=
> "*"/[f.factor^(f.exponent rem 2) for f in factorList smpsqfr])<br=
> >
> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 u : RR :=3D retrac=
> t(coerce(unit smpsqfr)@F) 
> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 sqfru :=3D squareF=
> ree u 
> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 ua :=3D "*&qu=
> ot;/[f.factor^(f.exponent quo 2) for f in factorList sqfru] 
> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 ub :=3D "*&qu=
> ot;/[f.factor^(f.exponent rem 2) for f in factorList sqfru] 
> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 ua*a*sqrt(ub*b)<br=
> >
> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 K =3D=3D> Kerne=
> l F 
> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 PP =3D=3D> Spar=
> seMultivariatePolynomial(RR, K) 
> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 res :=3D froot(s, =
> 2)$PolynomialRoots(IndexedExponents K, K, RR, PP, F) 
> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 res.coef * (res.ra=
> dicand)^(1/res.exponent) 
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0else 
> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0s^(1/2) </=
> div>
>
> 
>
> -- 
> You received this message because you are subscribed to the Google Groups &=
> quot;FriCAS - computer algebra system" group. 
> To unsubscribe from this group and stop receiving emails from it, send an e=
> mail to <a href=3D"mailto:fricas-devel...@googlegroups.com">fricas=
> -devel+un...@googlegroups.com</a>. 
> To post to this group, send email to <a href=3D"mailto:fricas-devel@googleg=
> roups.com">fricas...@googlegroups.com</a>. 
> To view this discussion on the web visit <a href=3D"https://groups.google.c=
> om/d/msgid/fricas-devel/CAGBJN93gf%3Du2J8MkDi3QPvNC7wVn2BHAfEnXUcM6-UHM_7DZ=
> WA%40mail.gmail.com?utm_medium=3Demail&utm_source=3Dfooter">https://groups.=
> google.com/d/msgid/fricas-devel/CAGBJN93gf%3Du2J8MkDi3QPvNC7wVn2BHAfEnXUcM6=
> -UHM_7DZWA%40mail.gmail.com</a>. 
> For more options, visit <a href=3D"https://groups.google.com/d/optout">http=
> s://groups.google.com/d/optout</a>. 
>
> --0000000000004e2b0f05867325e0--
>

--
Waldek Hebisch

Reply all

Reply to author

Forward