| blob | fa88889815ae18235962c6680846a86ddfb3ec1b |
1 /* Original version of this file copyright 1999 by Michael Wester,
2 * and retrieved from http://www.math.unm.edu/~wester/demos/Algebra/problems.macsyma
3 * circa 2006-10-23.
4 *
5 * Released under the terms of the GNU General Public License, version 2,
6 * per message dated 2007-06-03 from Michael Wester to Robert Dodier
7 * (contained in the file wester-gpl-permission-message.txt).
8 *
9 * See: "A Critique of the Mathematical Abilities of CA Systems"
10 * by Michael Wester, pp 25--60 in
11 * "Computer Algebra Systems: A Practical Guide", edited by Michael J. Wester
12 * and published by John Wiley and Sons, Chichester, United Kingdom, 1999.
13 */
14 /* ----------[ M a c s y m a ]---------- */
15 /* ---------- Initialization ---------- */
16 showtime: all$
17 prederror: false$
18 /* ---------- Algebra ---------- */
19 /* One would think that the simplification 2 2^n => 2^(n + 1) would happen
20 automatically or at least easily ... */
21 2*2^n;
22 subst(a = 2, subst(2 = a, %));
23 /* And how about 4 2^n => 2^(n + 2)? [Richard Fateman] */
24 4*2^n;
25 map(factor, %);
26 /* (-1)^[n (n + 1)] => 1 for integer n */
27 declare(n, integer)$
28 (-1)^(n*(n + 1));
29 remove(n, integer)$
30 /* Also easy => 2 (3 x - 5) */
31 factor(6*x - 10);
32 /* Univariate gcd: gcd(p1, p2) => 1, gcd(p1 q, p2 q) => q [Richard Liska] */
33 p1: 64*x^34 - 21*x^47 - 126*x^8 - 46*x^5 - 16*x^60 - 81$
34 p2: 72*x^60 - 25*x^25 - 19*x^23 - 22*x^39 - 83*x^52 + 54*x^10 + 81$
35 q: 34*x^19 - 25*x^16 + 70*x^7 + 20*x^3 - 91*x - 86$
36 gcd(p1, p2);
37 gcd(expand(p1*q), expand(p2*q)) - q;
38 /* resultant(p1 q, p2 q) => 0 */
39 resultant(expand(p1*q), expand(p2*q), x);
40 /* How about factorization? => p1 * p2 */
41 factor(expand(p1 * p2));
42 remvalue(p1, p2, q)$
43 /* Multivariate gcd: gcd(p1, p2) => 1, gcd(p1 q, p2 q) => q */
44 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$
45 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$
46 q: 11*x^12*y^7*z^13 - 23*x^2*y^8*z^10 + 47*x^17*y^5*z^8$
47 gcd(p1, p2);
48 gcd(expand(p1*q), expand(p2*q)) - q;
49 /* How about factorization? => p1 * p2 */
50 factor(expand(p1 * p2));
51 remvalue(p1, p2, q)$
52 /* => x^n for n > 0 [Chris Hurlburt] */
53 gcd(2*x^(n + 4) - x^(n + 2), 4*x^(n + 1) + 3*x^n);
54 /* Resultants. If the resultant of two polynomials is zero, this implies they
55 have a common factor. See Keith O. Geddes, Stephen R. Czapor and George
56 Labahn, _Algorithms for Computer Algebra_, Kluwer Academic Publishers, 1992,
57 p. 286 => 0 */
58 resultant(3*x^4 + 3*x^3 + x^2 - x - 2, x^3 - 3*x^2 + x + 5, x);
59 /* Numbers are nice, but symbols allow for variability---try some high school
60 algebra: rational simplification => (x - 2)/(x + 2) */
61 (x^2 - 4)/(x^2 + 4*x + 4);
62 ratsimp(%);
63 /* This example requires more sophistication => e^(x/2) - 1 */
64 [(%e^x - 1)/(%e^(x/2) + 1), (exp(x) - 1)/(exp(x/2) + 1)];
65 radcan(%);
66 /* Expand and factor polynomials */
67 (x + 1)^20;
68 expand(%);
69 diff(%, x);
70 factor(%);
71 /* Completely factor this polynomial, then try to multiply it back together! */
72 solve(x^3 + x^2 - 7 = 0, x);
73 apply("*", map(lambda([e], lhs(e) - rhs(e)), %));
74 ratsimp(expand(%));
75 x^100 - 1;
76 factor(%);
77 /* Factorization over the complex rationals
78 => (2 x + 3 i) (2 x - 3 i) (x + 1 + 4 i) (x + 1 - 4 i) */
79 gfactor(4*x^4 + 8*x^3 + 77*x^2 + 18*x + 153);
80 /* Algebraic extensions */
81 algebraic: true$
82 tellrat(sqrt2^2 - 2);
83 /* => sqrt2 + 1 */
84 rat(1/(sqrt2 - 1));
85 /* => (x^2 - 2 x - 3)/(x - sqrt2) = (x + 1) (x - 3)/(x - sqrt2)
86 [Richard Liska] */
87 (x^3 + (sqrt2 - 2)*x^2 - (2*sqrt2 + 3)*x - 3*sqrt2)/(x^2 - 2);
88 rat(%);
89 factor(%);
90 factor(%, sqrt2^2 - 2);
91 untellrat(sqrt2)$
92 /* Multiple algebraic extensions */
93 tellrat(sqrt3^2 - 3, cbrt2^3 - 2);
94 /* => 2 cbrt2 + 8 sqrt3 + 18 cbrt2^2 + 12 cbrt2 sqrt3 + 9 */
95 rat((cbrt2 + sqrt3)^4);
96 untellrat(sqrt3, cbrt2)$
97 algebraic: false$
98 /* Factor polynomials over finite fields and field extensions */
99 p: x^4 - 3*x^2 + 1;
100 factor(p);
101 /* => (x - 2)^2 (x + 2)^2 mod 5 */
102 ev(factor(p), modulus:5);
103 expand(%);
104 /* => (x^2 + x + 1) (x^9 - x^8 + x^6 - x^5 + x^3 - x^2 + 1) mod 65537
105 [Paul Zimmermann] */
106 ev(factor(x^11 + x + 1), modulus:65537);
107 /* => (x - phi) (x + phi) (x - phi + 1) (x + phi - 1)
108 where phi^2 - phi - 1 = 0 or phi = (1 +- sqrt(5))/2 */
109 factor(p, phi^2 - phi - 1);
110 remvalue(p)$
111 expand((x - 2*y^2 + 3*z^3)^20)$
112 factor(%);
113 expand((sin(x) - 2*cos(y)^2 + 3*tan(z)^3)^20)$
114 factor(%);
115 /* expand[(1 - c^2)^5 (1 - s^2)^5 (c^2 + s^2)^10] => c^10 s^10 when
116 c^2 + s^2 = 1 [modification of a problem due to Richard Liska] */
117 expand((1 - c^2)^5 * (1 - s^2)^5 * (c^2 + s^2)^10)$
118 grobner([%, c^2 + s^2 - 1]);
119 factor(%);
120 /* => (x + y) (x - y) mod 3 */
121 ev(factor(4*x^2 - 21*x*y + 20*y^2), modulus:3);
122 /* => 1/4 (x + y) (2 x + y [-1 + i sqrt(3)]) (2 x + y [-1 - i sqrt(3)]) */
123 factor(x^3 + y^3, isqrt3^2 + 3);
124 /* Partial fraction decomposition => 3/(x + 2) - 2/(x + 1) + 2/(x + 1)^2 */
125 (x^2 + 2*x + 3)/(x^3 + 4*x^2 + 5*x + 2);
126 partfrac(%, x);
127 /* Noncommutative algebra: note that (A B C)^(-1) = C^(-1) B^(-1) A^(-1)
128 => A B C A C B - C^(-1) B^(-1) C B */
129 (A.B.C - (A.B.C)^^(-1)) . A.C.B;
130 expand(%);
131 /* Jacobi's identity: [A, B, C] + [B, C, A] + [C, A, B] = 0 where [A, B, C] =
132 [A, [B, C]] and [A, B] = A B - B A is the commutator of A and B */
133 comm2(A, B):= A . B - B . A$
134 comm3(A, B, C):= comm2(A, comm2(B, C))$
135 comm2(A, B);
136 comm3(A, B, C) + comm3(B, C, A) + comm3(C, A, B);
137 expand(%);
138 remfunction(comm2, comm3)$
139 /* ---------- Quit ---------- */
140 quit();