1% polydiv.tst -*- REDUCE -*- 2 3% Test and demonstration file for enhanced polynomial division 4% file polydiv.red. 5 6% F.J.Wright@Maths.QMW.ac.uk, 7 Nov 1995. 7 8% The example from "Computer Algebra" by Davenport, Siret & Tournier, 9% first edition, section 2.3.3. 10 11% First check that remainder still works as before. 12 13% Compute the gcd of the polynomials a and b by Euclid's algorithm: 14a := aa := x^8 + x^6 - 3x^4 - 3x^3 + 8x^2 + 2x - 5; 15b := bb := 3x^6 + 5x^4 - 4x^2 - 9x + 21; 16on rational; off allfac; 17c := remainder(a, b); a := b$ b := c$ 18c := remainder(a, b); a := b$ b := c$ 19c := remainder(a, b); a := b$ b := c$ 20c := remainder(a, b); a := b$ b := c$ 21c := remainder(a, b); 22off rational; 23 24% Repeat using pseudo-remainders, to avoid rational arithmetic: 25a := aa; 26b := bb; 27c := pseudo_remainder(a, b); a := b$ b := c$ 28c := pseudo_remainder(a, b); a := b$ b := c$ 29c := pseudo_remainder(a, b); a := b$ b := c$ 30c := pseudo_remainder(a, b); a := b$ b := c$ 31c := pseudo_remainder(a, b); 32 33 34% Example from Chris Herssens <herc@sulu.luc.ac.be> 35% involving algebraic numbers in the coefficient ring 36% (for which naive pseudo-division fails in REDUCE): 37factor x; 38a:=8*(15*sqrt(2)*x**3 + 18*sqrt(2)*x**2 + 10*sqrt(2)*x + 12*sqrt(2) - 39 5*x**4 - 6*x**3 - 30*x**2 - 36*x); 40b:= - 16320*sqrt(2)*x**3 - 45801*sqrt(2)*x**2 - 50670*sqrt(2)*x - 41 26534*sqrt(2) + 15892*x**3 + 70920*x**2 + 86352*x + 24780; 42 43pseudo_remainder(a, b, x); 44% Note: We must specify the division variable even though the 45% polynomials are apparently univariate: 46pseudo_remainder(a, b); 47 48% Confirm that quotient * b + remainder = constant * a: 49pseudo_divide(a, b, x); 50first ws * b + second ws; 51ws / a; % is this constant? 52on rationalize; 53ws; % yes, it is constant 54off rationalize; 55 56on allfac; remfac x; 57 58procedure test_pseudo_division(a, b, x); 59 begin scalar qr, L; 60 qr := pseudo_divide(a, b, x); 61 L := lcof(b,x); 62 %% For versions of REDUCE prior to 3.6 use: 63 %% L := if b freeof x then b else lcof(b,x); 64 if first qr * b + second qr = 65 L^(deg(a,x)-deg(b,x)+1) * a then 66 write "Pseudo-division OK" 67 else 68 write "Pseudo-division failed" 69 end; 70 71a := 5x^4 + 4x^3 + 3x^2 + 2x + 1; 72test_pseudo_division(a, x, x); 73test_pseudo_division(a, x^3, x); 74test_pseudo_division(a, x^5, x); 75test_pseudo_division(a, x^3 + x, x); 76test_pseudo_division(a, 0, x); % intentional error! 77test_pseudo_division(a, 1, x); 78 79test_pseudo_division(5x^3 + 7y^2, 2x - y, x); 80test_pseudo_division(5x^3 + 7y^2, 2x - y, y); 81 82end; 83