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; 15 16 17 8 6 4 3 2 18a := aa := x + x - 3*x - 3*x + 8*x + 2*x - 5 19 20b := bb := 3x^6 + 5x^4 - 4x^2 - 9x + 21; 21 22 23 6 4 2 24b := bb := 3*x + 5*x - 4*x - 9*x + 21 25 26on rational; 27 28 off allfac; 29 30 31c := remainder(a, b); 32 33 34 5 4 1 2 1 35c := - ---*x + ---*x - --- 36 9 9 3 37 a := b$ 38 39 b := c$ 40 41 42c := remainder(a, b); 43 44 45 117 2 441 46c := - -----*x - 9*x + ----- 47 25 25 48 a := b$ 49 50 b := c$ 51 52 53c := remainder(a, b); 54 55 56 233150 102500 57c := --------*x - -------- 58 19773 6591 59 a := b$ 60 61 b := c$ 62 63 64c := remainder(a, b); 65 66 67 1288744821 68c := - ------------ 69 543589225 70 a := b$ 71 72 b := c$ 73 74 75c := remainder(a, b); 76 77 78c := 0 79 80off rational; 81 82 83 84% Repeat using pseudo-remainders, to avoid rational arithmetic: 85a := aa; 86 87 88 8 6 4 3 2 89a := x + x - 3*x - 3*x + 8*x + 2*x - 5 90 91b := bb; 92 93 94 6 4 2 95b := 3*x + 5*x - 4*x - 9*x + 21 96 97c := pseudo_remainder(a, b); 98 99 100 4 2 101c := - 15*x + 3*x - 9 102 a := b$ 103 104 b := c$ 105 106 107c := pseudo_remainder(a, b); 108 109 110 2 111c := 15795*x + 30375*x - 59535 112 a := b$ 113 114 b := c$ 115 116 117c := pseudo_remainder(a, b); 118 119 120c := 1254542875143750*x - 1654608338437500 121 a := b$ 122 123 b := c$ 124 125 126c := pseudo_remainder(a, b); 127 128 129c := 12593338795500743100931141992187500 130 a := b$ 131 132 b := c$ 133 134 135c := pseudo_remainder(a, b); 136 137 138c := 0 139 140 141 142% Example from Chris Herssens <herc@sulu.luc.ac.be> 143% involving algebraic numbers in the coefficient ring 144% (for which naive pseudo-division fails in REDUCE): 145factor x; 146 147 148a:=8*(15*sqrt(2)*x**3 + 18*sqrt(2)*x**2 + 10*sqrt(2)*x + 12*sqrt(2) - 149 5*x**4 - 6*x**3 - 30*x**2 - 36*x); 150 151 152 4 3 2 153a := - 40*x + x *(120*sqrt(2) - 48) + x *(144*sqrt(2) - 240) 154 155 + x*(80*sqrt(2) - 288) + 96*sqrt(2) 156 157b:= - 16320*sqrt(2)*x**3 - 45801*sqrt(2)*x**2 - 50670*sqrt(2)*x - 158 26534*sqrt(2) + 15892*x**3 + 70920*x**2 + 86352*x + 24780; 159 160 161 3 2 162b := x *( - 16320*sqrt(2) + 15892) + x *( - 45801*sqrt(2) + 70920) 163 164 + x*( - 50670*sqrt(2) + 86352) - 26534*sqrt(2) + 24780 165 166 167pseudo_remainder(a, b, x); 168 169 170 2 3/2 171x *( - 51343372800*2 + 72663731640*2 + 106394745600*sqrt(2) - 152808065280) + 172 173 3/2 174 x*( - 77924736000*2 + 111722451600*2 + 167518488000*sqrt(2) - 236076547200) 175 176 3/2 177 - 26395315200*2 + 21508247760*2 + 58006274400*sqrt(2) - 51393323520 178 179% Note: We must specify the division variable even though the 180% polynomials are apparently univariate: 181pseudo_remainder(a, b); 182 183 184*** Main division variable selected is 2**(1/2) 185 186 7 6 5 4 3 2 187652800*x + 708360*x - 2656800*x - 2660160*x + 4017600*x + 3676320*x 188 189 - 2630400*x - 2378880 190 191 192% Confirm that quotient * b + remainder = constant * a: 193pseudo_divide(a, b, x); 194 195 196{x*(652800*sqrt(2) - 635680) - 1958400*2 + 858360*sqrt(2) + 2073984, 197 198 2 3/2 199 x *( - 51343372800*2 + 72663731640*2 + 106394745600*sqrt(2) - 152808065280) 200 201 + x 202 203 3/2 204 *( - 77924736000*2 + 111722451600*2 + 167518488000*sqrt(2) - 236076547200) 205 206 3/2 207 - 26395315200*2 + 21508247760*2 + 58006274400*sqrt(2) - 51393323520} 208 209first ws * b + second ws; 210 211 212 4 213x *(20748595200*sqrt(2) - 31409618560) 214 215 3 216 + x *(119127169920*sqrt(2) - 162183113472) 217 218 2 219 + x *(237566198016*sqrt(2) - 337847596800) 220 221 + x*(212209122560*sqrt(2) - 309143634432) + 75383084544*sqrt(2) - 99593256960 222 223ws / a; 224 225 226 4 3 227(x *(2593574400*sqrt(2) - 3926202320) + x *(14890896240*sqrt(2) - 20272889184) 228 229 2 230 + x *(29695774752*sqrt(2) - 42230949600) 231 232 + x*(26526140320*sqrt(2) - 38642954304) + 9422885568*sqrt(2) - 12449157120)/( 233 234 4 3 2 235 - 5*x + x *(15*sqrt(2) - 6) + x *(18*sqrt(2) - 30) + x*(10*sqrt(2) - 36) 236 237 + 12*sqrt(2)) 238 % is this constant? 239on rationalize; 240 241 242ws; 243 244 245 - 518714880*sqrt(2) + 785240464 246 % yes, it is constant 247off rationalize; 248 249 250 251on allfac; 252 253 remfac x; 254 255 256 257procedure test_pseudo_division(a, b, x); 258 begin scalar qr, L; 259 qr := pseudo_divide(a, b, x); 260 L := lcof(b,x); 261 %% For versions of REDUCE prior to 3.6 use: 262 %% L := if b freeof x then b else lcof(b,x); 263 if first qr * b + second qr = 264 L^(deg(a,x)-deg(b,x)+1) * a then 265 write "Pseudo-division OK" 266 else 267 write "Pseudo-division failed" 268 end; 269 270 271test_pseudo_division 272 273 274a := 5x^4 + 4x^3 + 3x^2 + 2x + 1; 275 276 277 4 3 2 278a := 5*x + 4*x + 3*x + 2*x + 1 279 280test_pseudo_division(a, x, x); 281 282 283Pseudo-division OK 284 285test_pseudo_division(a, x^3, x); 286 287 288Pseudo-division OK 289 290test_pseudo_division(a, x^5, x); 291 292 293Pseudo-division OK 294 295test_pseudo_division(a, x^3 + x, x); 296 297 298Pseudo-division OK 299 300test_pseudo_division(a, 0, x); 301 302 303***** Zero divisor 304 % intentional error! 305test_pseudo_division(a, 1, x); 306 307 308Pseudo-division OK 309 310 311test_pseudo_division(5x^3 + 7y^2, 2x - y, x); 312 313 314Pseudo-division OK 315 316test_pseudo_division(5x^3 + 7y^2, 2x - y, y); 317 318 319Pseudo-division OK 320 321 322end; 323 324Tested on x86_64-pc-windows CSL 325Time (counter 1): 0 ms 326 327End of Lisp run after 0.00+0.04 seconds 328real 0.20 329user 0.00 330sys 0.06 331