Mon Feb 23 21:58:07 MST 1998 aquarius% mupad *----* MuPAD 1.4.0 -- Multi Processing Algebra Data Tool /| /| *----* | Copyright (c) 1997 - 98 by SciFace Software GmbH | *--|-* All rights reserved. |/ |/ *----* Licensed to: Michael Wester >> # ----------[ M u P A D ]---------- # >> # ---------- Initialization ---------- # >> TEXTWIDTH:= 80: >> read("../../Time.mupad"): >> # ---------- Algebra ---------- # >> # One would think that the simplification 2 2^n => 2^(n + 1) would happen &> automatically or at least easily ... # >> 2*2^n; n 2 2 Time: 90 msec Type: "_mult" >> combine(%); n + 1 2 Time: 580 msec Type: "_power" >> # And how about 4 2^n => 2^(n + 2)? [Richard Fateman] # >> 4*2^n; n 4 2 Time: 70 msec Type: "_mult" >> # (-1)^[n (n + 1)] => 1 for integer n # >> assume(n, Type::Integer): >> (-1)^(n*(n + 1)); n (n + 1) (-1) Time: 6630 msec Type: "_power" >> simplify(%); n (n + 1) (-1) Time: 990 msec Type: "_power" >> unassume(n, Type::Integer): >> # Also easy => 2 (3 x - 5) # >> Factor(6*x - 10); 6 x - 10 Time: 880 msec Type: "_plus" >> factor(6*x - 10); [2, 3 x - 5, 1] Time: 130 msec Type: DOM_LIST >> # Univariate gcd: gcd(p1, p2) => 1, gcd(p1 q, p2 q) => q [Richard Liska] # >> p1:= 64*x^34 - 21*x^47 - 126*x^8 - 46*x^5 - 16*x^60 - 81: >> p2:= 72*x^60 - 25*x^25 - 19*x^23 - 22*x^39 - 83*x^52 + 54*x^10 + 81: >> q:= 34*x^19 - 25*x^16 + 70*x^7 + 20*x^3 - 91*x - 86: >> gcd(p1, p2); 1 Time: 1020 msec Type: DOM_INT >> gcd(expand(p1*q), expand(p2*q)) - q; 0 Time: 230 msec Type: DOM_INT >> # resultant(p1 q, p2 q) => 0 # >> resultant(expand(p1*q), expand(p2*q), x); 0 Time: 7620 msec Type: DOM_INT >> # How about factorization? => p1 * p2 # >> Factor(expand(p1 * p2)); 5 8 34 47 60 - (46 x + 126 x - 64 x + 21 x + 16 x + 81) 10 23 25 39 52 60 (54 x - 19 x - 25 x - 22 x - 83 x + 72 x + 81) Time: 89020 msec Type: "_mult" >> p1:= NIL: >> p2:= NIL: >> q:= NIL: >> # Multivariate gcd: gcd(p1, p2) => 1, gcd(p1 q, p2 q) => q # >> p1:= 24*x*y^19*z^8 - 47*x^17*y^5*z^8 + 6*x^15*y^9*z^2 - 3*x^22 + 5: >> p2:= 34*x^5*y^8*z^13 + 20*x^7*y^7*z^7 + 12*x^9*y^16*z^4 + 80*y^14*z: >> q:= 11*x^12*y^7*z^13 - 23*x^2*y^8*z^10 + 47*x^17*y^5*z^8: >> gcd(p1, p2); 1 Time: 1450 msec Type: DOM_INT >> gcd(expand(p1*q), expand(p2*q)) - q; 0 Time: 2590 msec Type: DOM_INT >> # How about factorization? => p1 * p2 # >> Factor(expand(p1 * p2)); 7 7 7 6 5 12 9 9 3 2 y z (40 y + 10 x z + 17 x y z + 6 x y z ) 19 8 22 15 9 2 17 5 8 (24 x y z - 3 x + 6 x y z - 47 x y z + 5) Time: 30680 msec Type: "_mult" >> p1:= NIL: >> p2:= NIL: >> q:= NIL: >> # => x^n for n > 0 [Chris Hurlburt] # >> assume(n > 0): >> gcd(2*x^(n + 4) - x^(n + 2), 4*x^(n + 1) + 3*x^n); 1 Time: 6900 msec Type: DOM_INT >> unassume(n > 0): >> # Resultants. If the resultant of two polynomials is zero, this implies they &> have a common factor. See Keith O. Geddes, Stephen R. Czapor and George &> Labahn, _Algorithms for Computer Algebra_, Kluwer Academic Publishers, 1992, &> p. 286 => 0 # >> resultant(3*x^4 + 3*x^3 + x^2 - x - 2, x^3 - 3*x^2 + x + 5, x); 0 Time: 520 msec Type: DOM_INT >> # Numbers are nice, but symbols allow for variability---try some high school &> algebra: rational simplification => (x - 2)/(x + 2) # >> (x^2 - 4)/(x^2 + 4*x + 4); 2 x - 4 ------------ 2 4 x + x + 4 Time: 230 msec Type: "_mult" >> normal(%); x - 2 ----- x + 2 Time: 280 msec Type: "_mult" >> # This example requires more sophistication => e^(x/2) - 1 # >> [(E^x - 1)/(E^(x/2) + 1), (exp(x) - 1)/(exp(x/2) + 1)]; -- x -- | exp(1) - 1 exp(x) - 1 | | -----------, ------------ | | x / x \ | | - exp| - | + 1 | | 2 \ 2 / | -- exp(1) + 1 -- Time: 240 msec Type: DOM_LIST >> simplify(%, exp); -- x -- | exp(1) - 1 / x \ | | ----------------, exp| - | - 1 | | x 1/2 \ 2 / | -- (exp(1) ) + 1 -- Time: 410 msec Type: DOM_LIST >> # Expand and factor polynomials # >> (x + 1)^20; 20 (x + 1) Time: 230 msec Type: "_power" >> expand(%); 2 3 4 5 6 7 8 20 x + 190 x + 1140 x + 4845 x + 15504 x + 38760 x + 77520 x + 125970 x + 9 10 11 12 13 14 167960 x + 184756 x + 167960 x + 125970 x + 77520 x + 38760 x + 15 16 17 18 19 20 15504 x + 4845 x + 1140 x + 190 x + 20 x + x + 1 Time: 240 msec Type: "_plus" >> diff(%, x); 2 3 4 5 6 7 380 x + 3420 x + 19380 x + 77520 x + 232560 x + 542640 x + 1007760 x + 8 9 10 11 12 1511640 x + 1847560 x + 1847560 x + 1511640 x + 1007760 x + 13 14 15 16 17 18 19 542640 x + 232560 x + 77520 x + 19380 x + 3420 x + 380 x + 20 x + 20 Time: 240 msec Type: "_plus" >> Factor(%); 19 20 (x + 1) Time: 470 msec Type: "_mult" >> # Completely factor this polynomial, then try to multiply it back together! # >> solve(x^3 + x^2 - 7 = 0, x); { / 1/2 1/2 \1/3 / 1/2 1/2 \1/3 { | 108 1295 | | 108 1295 | { | 187/54 - -------------- | | -------------- + 187/54 | { \ 108 / \ 108 / { - ------------------------------ - ------------------------------ - { 2 2 / / 1/2 1/2 \1/3 / 1/2 1/2 \1/3 1/2 | | 108 1295 | | 108 1295 | 1/2 I 3 | | -------------- + 187/54 | - | 187/54 - -------------- | \ \ 108 / \ 108 / \ / / 1/2 1/2 \1/3 | 1/2 | | 108 1295 | | - 1/3, 1/2 I 3 | | -------------- + 187/54 | - / \ \ 108 / / 1/2 1/2 \1/3 | 108 1295 | / 1/2 1/2 \1/3 \ | 187/54 - -------------- | | 108 1295 | | \ 108 / | 187/54 - -------------- | | - ------------------------------ - \ 108 / / 2 / 1/2 1/2 \1/3 | 108 1295 | | -------------- + 187/54 | / 1/2 1/2 \1/3 \ 108 / | 108 1295 | ------------------------------ - 1/3, | 187/54 - -------------- | + 2 \ 108 / } } / 1/2 1/2 \1/3 } | 108 1295 | } | -------------- + 187/54 | - 1/3 } \ 108 / } Time: 2380 msec Type: DOM_SET >> _mult(op(map(%, func(x - e, e)))); / | / / 1/2 1/2 \1/3 / 1/2 1/2 \1/3 \ | | | 108 1295 | | 108 1295 | | | | x - | 187/54 - -------------- | - | -------------- + 187/54 | + 1/3 | | \ \ 108 / \ 108 / / \ / 1/2 1/2 \1/3 / 1/2 1/2 \1/3 | 108 1295 | | 108 1295 | | 187/54 - -------------- | | -------------- + 187/54 | \ 108 / \ 108 / x + ------------------------------ + ------------------------------ - 1/2 I 2 2 / / 1/2 1/2 \1/3 / 1/2 1/2 \1/3 \ 1/2 | | 108 1295 | | 108 1295 | | 3 | | -------------- + 187/54 | - | 187/54 - -------------- | | + \ \ 108 / \ 108 / / \ / / 1/2 1/2 \1/3 / 1/2 1/2 \1/3 | | | 108 1295 | | 108 1295 | | | | 187/54 - -------------- | | -------------- + 187/54 | | | \ 108 / \ 108 / 1/3 | | x + ------------------------------ + ------------------------------ + / \ 2 2 / / 1/2 1/2 \1/3 / 1/2 1/2 \1/3 1/2 | | 108 1295 | | 108 1295 | 1/2 I 3 | | -------------- + 187/54 | - | 187/54 - -------------- | \ \ 108 / \ 108 / \ | \ | | | | + 1/3 | / / Time: 400 msec Type: "_mult" >> simplify(expand(%)); / / 1/2 \ / 1/2 \ \1/3 x 2 3 | | 3885 | | 3885 | | - + x + x - | | 187/54 - ------- | | ------- + 187/54 | | - 3 \ \ 18 / \ 18 / / / / 1/2 \ / 1/2 \ \1/3 | | 3885 | | 3885 | | 3 x | | 187/54 - ------- | | ------- + 187/54 | | - 62/9 \ \ 18 / \ 18 / / Time: 2640 msec Type: "_plus" >> float(%); 2 3 0.000000000000000006125742266 x + x + x - 7.0 Time: 280 msec Type: "_plus" >> x^100 - 1; 100 x - 1 Time: 280 msec Type: "_plus" >> Factor(%); 2 2 3 4 5 10 15 20 (x - 1) (x + 1) (x + 1) (x + x + x + x + 1) (x + x + x + x + 1) 2 3 4 4 2 6 8 10 5 15 20 (x - x - x + x + 1) (x - x - x + x + 1) (x - x - x + x + 1) 20 10 30 40 (x - x - x + x + 1) Time: 410 msec Type: "_mult" >> # Factorization over the complex rationals &> => (2 x + 3 i) (2 x - 3 i) (x + 1 + 4 i) (x + 1 - 4 i) # >> QI:= Dom::AlgebraicExtension(Dom::Rational, i^2 + 1); 2 Dom::AlgebraicExtension(Dom::Rational, i + 1 = 0, i) Time: 1700 msec Type: DOM_DOMAIN >> QI::name:= "QI": >> Factor(poly(4*x^4 + 8*x^3 + 77*x^2 + 18*x + 153, QI)); / 3 i \ / 3 i \ 4 poly| x + ---, [x], QI | poly(x + (4 i + 1), [x], QI) poly| x - ---, [x], QI | \ 2 / \ 2 / poly(x + (1 - 4 i), [x], QI) Time: 3580 msec Type: "_mult" >> # Algebraic extensions # >> Qsqrt2:= Dom::AlgebraicExtension(Dom::Rational, sqrt2^2 - 2): >> Qsqrt2::name:= "Qsqrt2": >> # => sqrt2 + 1 # >> 1/(sqrt2 - 1); 1 --------- sqrt2 - 1 Time: 920 msec Type: "_power" >> # => (x^2 - 2 x - 3)/(x - sqrt2) = (x + 1) (x - 3)/(x - sqrt2) &> [Richard Liska] # >> p:= (x^3 + (sqrt2 - 2)*x^2 - (2*sqrt2 + 3)*x - 3*sqrt2)/(x^2 - 2); 3 2 x - 3 sqrt2 - x (2 sqrt2 + 3) + x (sqrt2 - 2) ----------------------------------------------- 2 x - 2 Time: 300 msec Type: "_mult" >> Factor(poly(p, Qsqrt2)); FAIL Time: 860 msec Type: DOM_FAIL >> [Factor(poly(numer(p), Qsqrt2)), Factor(poly(denom(p), Qsqrt2))]; [poly(x + 1, [sqrt2, x], Qsqrt2) poly(x - 3, [sqrt2, x], Qsqrt2) poly(sqrt2 + x, [sqrt2, x], Qsqrt2), poly(x - sqrt2, [x], Qsqrt2) poly(x + sqrt2, [x], Qsqrt2)] Time: 2820 msec Type: DOM_LIST >> %[1]/%[2]; Error: Argument out of range [_power] >> p:= NIL: >> sqrt2:= NIL: >> # Multiple algebraic extensions # >> Qsqrt3:= Dom::AlgebraicExtension(Dom::Rational, sqrt3^2 - 3): >> Qsqrt3::name:= "Qsqrt3": >> Qcbrt2:= Dom::AlgebraicExtension(Dom::Rational, cbrt2^3 - 2): >> Qcbrt2::name:= "Qcbrt2": >> # => 2 cbrt2 + 8 sqrt3 + 18 cbrt2^2 + 12 cbrt2 sqrt3 + 9 # >> (cbrt2 + sqrt3)^4; 4 (cbrt2 + sqrt3) Time: 2520 msec Type: "_power" >> expand(%); 4 4 3 3 2 2 cbrt2 + sqrt3 + 4 cbrt2 sqrt3 + 4 cbrt2 sqrt3 + 6 cbrt2 sqrt3 Time: 310 msec Type: "_plus" >> sqrt3:= NIL: >> cbrt2:= NIL: >> # Factor polynomials over finite fields and field extensions # >> p:= x^4 - 3*x^2 + 1; 4 2 x - 3 x + 1 Time: 900 msec Type: "_plus" >> Factor(p); 2 2 (x + x - 1) (x - x - 1) Time: 700 msec Type: "_mult" >> # => (x - 2)^2 (x + 2)^2 mod 5 # >> Factor(poly(p, IntMod(5))); 2 2 poly(x - 2, [x], IntMod(5)) poly(x + 2, [x], IntMod(5)) Time: 470 msec Type: "_mult" >> # => (x^2 + x + 1) (x^9 - x^8 + x^6 - x^5 + x^3 - x^2 + 1) mod 65537 &> [Paul Zimmermann] # >> Factor(poly(x^11 + x + 1, IntMod(65537))); 2 poly(x + x + 1, [x], IntMod(65537)) poly( 9 8 6 5 3 2 x + (-1) x + x + (-1) x + x + (-1) x + 1, [x], IntMod(65537)) Time: 600 msec Type: "_mult" >> # => (x - phi) (x + phi) (x - phi + 1) (x + phi - 1) &> where phi^2 - phi - 1 = 0 or phi = (1 +- sqrt(5))/2 # >> Qphi:= Dom::AlgebraicExtension(Dom::Rational, Phi^2 - Phi - 1): >> Qphi::name:= "Qphi": >> Factor(poly(p, Qphi)); poly(x + Phi, [x], Qphi) poly(x + (Phi - 1), [x], Qphi) poly(x - Phi, [x], Qphi) poly(x + (1 - Phi), [x], Qphi) Time: 3750 msec Type: "_mult" >> p:= NIL: >> expand((x - 2*y^2 + 3*z^3)^20): >> Factor(%); 2 3 20 (x - 2 y + 3 z ) Time: 6870 msec Type: "_power" >> expand((sin(x) - 2*cos(y)^2 + 3*tan(z)^3)^20): >> Factor(%); 3 3 2 3 20 (3 sin(z) + cos(z) sin(x) - 2 cos(y) cos(z) ) -------------------------------------------------- 60 cos(z) Time: 40640 msec Type: "_mult" >> # expand[(1 - c^2)^5 (1 - s^2)^5 (c^2 + s^2)^10] => c^10 s^10 when &> c^2 + s^2 = 1 [modification of a problem due to Richard Liska] # >> expand((1 - c^2)^5 * (1 - s^2)^5 * (c^2 + s^2)^10): >> groebner::gbasis([%, c^2 + s^2 - 1]); 12 10 14 16 18 20 2 2 [5 c - c - 10 c + 10 c - 5 c + c , c + s - 1] Time: 2240 msec Type: DOM_LIST >> map(%, Factor); 10 5 5 2 2 [c (c - 1) (c + 1) , c + s - 1] Time: 580 msec Type: DOM_LIST >> # => (x + y) (x - y) mod 3 # >> Factor(poly(4*x^2 - 21*x*y + 20*y^2, IntMod(3))); poly(x + (-1) y, [x, y], IntMod(3)) poly(x + y, [x, y], IntMod(3)) Time: 2770 msec Type: "_mult" >> # => 1/4 (x + y) (2 x + y [-1 + i sqrt(3)]) (2 x + y [-1 - i sqrt(3)]) # >> Qisqrt3:= Dom::AlgebraicExtension(Dom::Rational, isqrt3^2 + 3): >> Qisqrt3::name:= "Qisqrt3": >> Factor(poly(x^3 + y^3, Qisqrt3)); / / isqrt3 \ \ poly(x + y, [x, y], Qisqrt3) poly| x + | - ------ - 1/2 | y, [x, y], Qisqrt3 | \ \ 2 / / / / isqrt3 \ \ poly| x + | ------ - 1/2 | y, [x, y], Qisqrt3 | \ \ 2 / / Time: 4900 msec Type: "_mult" >> isqrt3:= NIL: >> # Partial fraction decomposition => 3/(x + 2) - 2/(x + 1) + 2/(x + 1)^2 # >> (x^2 + 2*x + 3)/(x^3 + 4*x^2 + 5*x + 2); 2 2 x + x + 3 ------------------- 2 3 5 x + 4 x + x + 2 Time: 740 msec Type: "_mult" >> partfrac(%); 2 2 3 -------- - ----- + ----- 2 x + 1 x + 2 (x + 1) Time: 600 msec Type: "_plus" >> # Noncommutative algebra: note that (A B C)^(-1) = C^(-1) B^(-1) A^(-1) &> => A B C A C B - C^(-1) B^(-1) C B # >> (A &ncmult B &ncmult C - (A &ncmult B &ncmult C)^(-1)) &> &ncmult A &ncmult C &ncmult B; / 1 \ ncmult| ncmult(A, B, C) - ---------------, A, C, B | \ ncmult(A, B, C) / Time: 390 msec Type: "function" >> # Jacobi's identity: [A, B, C] + [B, C, A] + [C, A, B] = 0 where [A, B, C] = &> [A, [B, C]] and [A, B] = A B - B A is the commutator of A and B # >> comm2:= proc(A, B) begin A &ncmult B - B &ncmult A end_proc: >> comm3:= proc(A, B, C) begin comm2(A, comm2(B, C)) end_proc: >> comm2(A, B); ncmult(A, B) - ncmult(B, A) Time: 1130 msec Type: "_plus" >> comm3(A, B, C) + comm3(B, C, A) + comm3(C, A, B); ncmult(C, ncmult(A, B) - ncmult(B, A)) - ncmult(ncmult(A, B) - ncmult(B, A), C) + ncmult(B, ncmult(C, A) - ncmult(A, C)) - ncmult(ncmult(C, A) - ncmult(A, C), B) + ncmult(A, ncmult(B, C) - ncmult(C, B)) - ncmult(ncmult(B, C) - ncmult(C, B), A) Time: 380 msec Type: "_plus" >> comm2:= NIL: >> comm3:= NIL: >> # ---------- Quit ---------- # >> quit real 243.56 user 241.99 sys 1.15