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