1% Tests of the poly package polynomial decomposition and gcds. 2 3 4% Test for the univariate and multivariate polynomial decomposition. 5 6% Herbert Melenk, ZIB Berlin, 1990. 7 8procedure testdecompose u; 9 begin scalar r,p,val,nextvar; 10 write "decomposition of ",u; 11 r := decompose u; 12 if length r = 1 then rederr "decomposition failed"; 13 write " leads to ",r; 14 % test if the result is algebraically correct. 15 r := reverse r; 16 nextvar := lhs first r; val := rhs first r; 17 r := rest r; 18 while not(r={}) do 19 << p := first r; r := rest r; 20 if 'equal = part(p,0) then 21 <<val := sub(nextvar=val,rhs p); nextvar := lhs p>> 22 else 23 val := sub(nextvar=val,p); 24 >>; 25 if val = u then write " O.K. " 26 else 27 <<write "**** reconstructed polynomial: "; 28 write val; 29 rederr "reconstruction leads to different polynomial"; 30 >>; 31 end; 32 33 34 % univariate decompositions 35testdecompose(x**4+x**2+1); 36testdecompose(x**6+9x**5+52x**4+177x**3+435x**2+630x+593); 37testdecompose(x**6+6x**4+x**3+9x**2+3x-5); 38testdecompose(x**8-88*x**7+2924*x**6-43912*x**5+263431*x**4-218900*x**3+ 39 65690*x**2-7700*x+234); 40 41 % multivariate cases 42testdecompose(u**2+v**2+2u*v+1); 43testdecompose(x**4+2x**3*y + 3x**2*y**2 + 2x*y**3 + y**4 + 2x**2*y 44 +2x*y**2 + 2y**3 + 5 x**2 + 5*x*y + 6*y**2 + 5y + 9); 45testdecompose sub(u=(2 x**2 + 17 x+y + y**3),u**2+2 u + 1); 46testdecompose sub(u=(2 x**2 *y + 17 x+y + y**3),u**2+2 u + 1); 47 48 % some cases which require a special (internal) mapping 49testdecompose ( (x + y)**2); 50testdecompose ((x + y**2)**2); 51testdecompose ( (x**2 + y)**2); 52testdecompose ( (u + v)**2 +10 ); 53 54 % the decomposition is not unique and might generate quite 55 % different images: 56testdecompose ( (u + v + 10)**2 -100 ); 57 58 % some special (difficult) cases 59testdecompose (X**4 + 88*X**3*Y + 2904*X**2*Y**2 - 10*X**2 60 + 42592*X*Y**3 - 440*X*Y + 234256*Y**4 - 4840*Y**2); 61 62 % a polynomial with complex coefficients 63on complex; 64testdecompose(X**4 + (88*I)*X**3*Y - 2904*X**2*Y**2 - 10*X**2 - 65 (42592*I)*X*Y**3 - (440*I)*X*Y + 234256*Y**4 + 4840*Y**2); 66off complex; 67 68 69% Examples given by J. Gutierrez and J.M. Olazabal. 70 71f1:=x**6-2x**5+x**4-3x**3+3x**2+5$ 72testdecompose(f1); 73 74f2:=x**32-1$ 75testdecompose(f2); 76 77f3:=x**4-(2/3)*x**3-(26/9)*x**2+x+3$ 78testdecompose(f3); 79 80f4:=sub(x=x**4-x**3-2x+1,x**3-x**2-1)$ 81testdecompose(f4); 82 83f5:=sub(x=f4,x**5-5)$ 84testdecompose(f5); 85 86clear f1,f2,f3,f4,f5; 87 88 89% Tests of gcd code. 90 91% The following examples were introduced in Moses, J. and Yun, D.Y.Y., 92% "The EZ GCD Algorithm", Proc. ACM 73 (1973) 159-166, and considered 93% further in Hearn, A.C., "Non-modular Computation of Polynomial GCD's 94% Using Trial Division", Proc. EUROSAM 79, 227-239, 72, published as 95% Lecture Notes on Comp. Science, # 72, Springer-Verlag, Berlin, 1979. 96 97on gcd; 98 99% The following is the best setting for this file. 100 101on ezgcd; 102 103% In systems that have the heugcd code, the following is also a 104% possibility, although not all examples complete in a reasonable time. 105 106% load heugcd; on heugcd; 107 108% The final alternative is to use neither ezgcd nor heugcd. In that case, 109% most examples take excessive amounts of computer time. 110 111share n; 112 113operator xx; 114 115% Case 1. 116 117for n := 2:5 118 do write gcd(((for i:=1:n sum xx(i))-1)*((for i:=1:n sum xx(i)) + 2), 119 ((for i:=1:n sum xx(i))+1) 120 *(-3xx(2)*xx(1)**2+xx(2)**2-1)**2); 121 122% Case 2. 123 124let d = (for i:=1:n sum xx(i)**n) + 1; 125 126for n := 2:7 do write gcd(d*((for i:=1:n sum xx(i)**n) - 2), 127 d*((for i:=1:n sum xx(i)**n) + 2)); 128 129 130for n := 2:7 do write gcd(d*((for i:=1:n sum xx(i)**n) - 2), 131 d*((for i:=1:n sum xx(i)**(n-1)) + 2)); 132 133% Case 3. 134 135let d = xx(2)**2*xx(1)**2 + (for i := 3:n sum xx(i)**2) + 1; 136 137for n := 2:5 138 do write gcd(d*(xx(2)*xx(1) + (for i:=3:n sum xx(i)) + 2)**2, 139 d*(xx(1)**2-xx(2)**2 + (for i:=3:n sum xx(i)**2) - 1)); 140 141% Case 4. 142 143let u = xx(1) - xx(2)*xx(3) + 1, 144 v = xx(1) - xx(2) + 3xx(3); 145 146gcd(u*v**2,v*u**2); 147 148gcd(u*v**3,v*u**3); 149 150gcd(u*v**4,v*u**4); 151 152gcd(u**2*v**4,v**2*u**4); 153 154 155% Case 5. 156 157let d = (for i := 1:n product (xx(i)+1)) - 3; 158 159for n := 2:5 do write gcd(d*for i := 1:n product (xx(i) - 2), 160 d*for i := 1:n product (xx(i) + 2)); 161 162clear d,u,v; 163 164 165% The following examples were discussed in Char, B.W., Geddes, K.O., 166% Gonnet, G.H., "GCDHEU: Heuristic Polynomial GCD Algorithm Based 167% on Integer GCD Computation", Proc. EUROSAM 84, 285-296, published as 168% Lecture Notes on Comp. Science, # 174, Springer-Verlag, Berlin, 1984. 169 170 171% Maple Problem 1. 172 173gcd(34*x**80-91*x**99+70*x**31-25*x**52+20*x**76-86*x**44-17*x**33 174 -6*x**89-56*x**54-17, 175 91*x**49+64*x**10-21*x**52-88*x**74-38*x**76-46*x**84-16*x**95 176 -81*x**72+96*x**25-20); 177 178% Maple Problem 2. 179 180g := 34*x**19-91*x+70*x**7-25*x**16+20*x**3-86; 181 182gcd(g * (64*x**34-21*x**47-126*x**8-46*x**5-16*x**60-81), 183 g * (72*x**60-25*x**25-19*x**23-22*x**39-83*x**52+54*x**10+81) ); 184 185% Maple Problem 3. 186 187gcd(3427088418+8032938293*x-9181159474*x**2-9955210536*x**3 188 +7049846077*x**4-3120124818*x**5-2517523455*x**6+5255435973*x**7 189 +2020369281*x**8-7604863368*x**9-8685841867*x**10+4432745169*x**11 190 -1746773680*x**12-3351440965*x**13-580100705*x**14+8923168914*x**15 191 -5660404998*x**16 +5441358149*x**17-1741572352*x**18 192 +9148191435*x**19-4940173788*x**20+6420433154*x**21+980100567*x**22 193 -2128455689*x**23+5266911072*x**24-8800333073*x**25-7425750422*x**26 194 -3801290114*x**27-7680051202*x**28-4652194273*x**29-8472655390*x**30 195 -1656540766*x**31+9577718075*x**32-8137446394*x**33+7232922578*x**34 196 +9601468396*x**35-2497427781*x**36-2047603127*x**37-1893414455*x**38 197 -2508354375*x**39-2231932228*x**40, 198 2503247071-8324774912*x+6797341645*x**2+5418887080*x**3 199 -6779305784*x**4+8113537696*x**5+2229288956*x**6+2732713505*x**7 200 +9659962054*x**8-1514449131*x**9+7981583323*x**10+3729868918*x**11 201 -2849544385*x**12-5246360984*x**13+2570821160*x**14-5533328063*x**15 202 -274185102*x**16+8312755945*x**17-2941669352*x**18-4320254985*x**19 203 +9331460166*x**20-2906491973*x**21-7780292310*x**22-4971715970*x**23 204 -6474871482*x**24-6832431522*x**25-5016229128*x**26-6422216875*x**27 205 -471583252*x**28+3073673916*x**29+2297139923*x**30+9034797416*x**31 206 +6247010865*x**32+5965858387*x**33-4612062748*x**34+5837579849*x**35 207 -2820832810*x**36-7450648226*x**37+2849150856*x**38+2109912954*x**39 208 +2914906138*x**40); 209 210% Maple Problem 4. 211 212g := 34271+80330*x-91812*x**2-99553*x**3+70499*x**4-31201*x**5 213 -25175*x**6+52555*x**7+20204*x**8-76049*x**9-86859*x**10; 214 215gcd(g * (44328-17468*x-33515*x**2-5801*x**3+89232*x**4-56604*x**5 216 +54414*x**6-17416*x**7+91482*x**8-49402*x**9+64205*x**10 217 +9801*x**11-21285*x**12+52669*x**13-88004*x**14-74258*x**15 218 -38013*x**16-76801*x**17-46522*x**18-84727*x**19-16565*x**20 219 +95778*x**21-81375*x**22+72330*x**23+96015*x**24-24974*x**25 220 -20476*x**26-18934*x**27-25084*x**28-22319*x**29+25033*x**30), 221 g * (-83248+67974*x+54189*x**2-67793*x**3+81136*x**4+22293*x**5 222 +27327*x**6+96600*x**7-15145*x**8+79816*x**9+37299*x**10 223 -28496*x**11-52464*x**12+25708*x**13-55334*x**14-2742*x**15 224 +83128*x**16-29417*x**17-43203*x**18+93315*x**19-29065*x**20 225 -77803*x**21-49717*x**22-64749*x**23-68325*x**24-50163*x**25 226 -64222*x**26-4716*x**27+30737*x**28+22972*x**29+90348*x**30)); 227 228% Maple Problem 5. 229 230gcd(-8472*x**4*y**10-8137*x**9*y**10-2497*x**4*y**4-2508*x**4*y**6 231 -8324*x**9*y**8-6779*x**9*y**6+2733*x**10*y**4+7981*x**7*y**3 232 -5246*x**6*y**2-274*x**10*y**3-4320, 233 15168*x**3*y-4971*x*y-2283*x*y**5+3074*x**6*y**10+6247*x**8*y**2 234 +2849*x**6*y**7-2039*x**7-2626*x**2*y**7+9229*x**6*y**5+2404*y**5 235 +1387*x**4*y**8+5602*x**5*y**2-6212*x**3*y**7-8561); 236 237% Maple Problem 6. 238 239g := -19*x**4*y**4+25*y**9+54*x*y**9+22*x**7*y**10-15*x**9*y**7-28; 240 241gcd(g*(91*x**2*y**9+10*x**4*y**8-88*x*y**3-76*x**2-16*x**10*y 242 +72*x**10*y**4-20), 243 g*(34*x**9-99*x**9*y**3-25*x**8*y**6-76*y**7-17*x**3*y**5 244 +89*x**2*y**8-17)); 245 246% Maple Problem 7. 247 248gcd(6713544209*x**9+8524923038*x**3*y**3*z**7+6010184640*x*z**7 249 +4126613160*x**3*y**4*z**9+2169797500*x**7*y**4*z**9 250 +2529913106*x**8*y**5*z**3+7633455535*y*z**3+1159974399*x**2*z**4 251 +9788859037*y**8*z**9+3751286109*x**3*y**4*z**3, 252 3884033886*x**6*z**8+7709443539*x*y**9*z**6 253 +6366356752*x**9*y**4*z**8+6864934459*x**3*y**2*z**6 254 +2233335968*x**4*y**9*z**3+2839872507*x**9*y**3*z 255 +2514142015*x*y*z**2+1788891562*x**4*y**6*z**6 256 +9517398707*x**8*y**7*z**2+7918789924*x**3*y*z**6 257 +6054956477*x**6*y**3*z**6); 258 259% Maple Problem 8. 260 261g := u**3*(x**2-y)*z**2+(u-3*u**2*x)*y*z-u**4*x*y+3; 262gcd(g * ((y**2+x)*z**2+u**5*(x*y+x**2)*z-y+5), 263 g * ((y**2-x)*z**2+u**5*(x*y-x**2)*z+y+9) ); 264 265% Maple Problem 9. 266 267g := 34*u**2*y**2*z-25*u**2*v*z**2-18*v*x**2*z**2-18*u**2*x**2*y*z+53 268 +x**3; 269gcd( g * (-85*u*v**2*y**2*z**2-25*u*v*x*y*z-84*u**2*v**2*y**2*z 270 +27*u**2*v*x**2*y**2*z-53*u*x*y**2*z+34*x**3), 271 g * (48*x**3-99*u*x**2*y**2*z-69*x*y*z-75*u*v*x*y*z**2 272 -43*u**2*v+91*u**2*v**2*y**2*z) ); 273 274% Maple Problem 10. 275 276gcd(-9955*v**9*x**3*y**4*z**8+2020*v*y**7*z**4 277 -3351*v**5*x**10*y**2*z**8-1741*v**10*x**2*y**9*z**6 278 -2128*v**8*y*z**3-7680*v**2*y**4*z**10-8137*v**9*x**10*y**4*z**4 279 -1893*v**4*x**4*y**6+6797*v**8*x*y**9*z**6 280 +2733*v**10*x**4*y**9*z**7-2849*v**2*x**6*y**2*z**5 281 +8312*v**3*x**3*y**10*z**3-7780*v**2*x*y*z**2 282 -6422*v**5*x**7*y**6*z**10+6247*v**8*x**2*y**8*z**3 283 -7450*v**7*x**6*y**7*z**4+3625*x**4*y**2*z**7+9229*v**6*x**5*y**6 284 -112*v**6*x**4*y**8*z**7-7867*v**5*x**8*y**5*z**2 285 -6212*v**3*x**7*z**5+8699*v**8*x**2*y**2*z**5 286 +4442*v**10*x**5*y**4*z+1965*v**10*y**3*z**3-8906*v**6*x*y**4*z**5 287 +5552*x**10*y**4+3055*v**5*x**3*y**6*z**2+6658*v**7*x**10*z**6 288 +3721*v**8*x**9*y**4*z**8+9511*v*x**6*y+5437*v**3*x**9*y**9*z**7 289 -1957*v**6*x**4*y*z**3+9214*v**3*x**9*y**3*z**7 290 +7273*v**2*x**8*y**4*z**10+1701*x**10*y**7*z**2 291 +4944*v**5*x**5*y**8*z**8-1935*v**3*x**6*y**10*z**7 292 +4029*x**6*y**10*z**3+9462*v**6*x**5*y**4*z**8-3633*v**4*x*y**7*z**5 293 -1876, 294 -5830*v**7*x**8*y*z**2-1217*v**8*x*y**2*z**5 295 -1510*v**9*x**3*y**10*z**10+7036*v**6*x**8*y**3*z**3 296 +1022*v**9*y**3*z**8+3791*v**8*x**3*y**7+6906*v**6*x*y*z**10 297 +117*v**7*x**2*y**4*z**4+6654*v**6*x**5*y**2*z**3 298 -7302*v**10*x**8*y**3-5343*v**8*x**5*y**9*z 299 -2244*v**9*x**3*y**8*z**9-3719*v**5*x**10*y**6*z**8 300 +2629*x**3*y**2*z**10+8517*x**9*y**6*z**7-9551*v**5*x**6*y**6*z**2 301 -7750*x**10*y**7*z**4-5035*v**5*x**2*y**5*z-5967*v**9*x**5*y**9*z**5 302 -8517*v**3*x**2*y**7*z**6-2668*v**10*y**9*z**4+1630*v**5*x**5*y*z**8 303 +9099*v**7*x**9*y**4*z**3-5358*v**9*x**5*y**6*z**2 304 +5766*v**5*y**3*z**4-3624*v*x**4*y**10*z**10 305 +8839*v**6*x**9*y**10*z**4+3378*x**7*y**2*z**5+7582*v**7*x*y**8*z**7 306 -85*v*x**2*y**9*z**6-9495*v**9*x**10*y**6*z**3+1983*v**9*x**3*y 307 -4613*v**10*x**4*y**7*z**6+5529*v**10*x*y**6 308 +5030*v**4*x**5*y**4*z**9-9202*x**6*y**3*z**9 309 -4988*v**2*x**2*y**10*z**4-8572*v**9*x**7*y**10*z**10 310 +4080*v**4*x**8*z**8-382*v**9*x**9*y**2*z**2-7326); 311 312end; 313 314