> # Worksheet to execute the examples in Kaltofen, May, Yang, & Zhi: > > # 'Approximate Factorization of Multivariate Polynomials Using > Singular Value Decomposition' > > restart; > > # Factorization routines: > read("multifac_1.3.mpl"): read("backwards_error.mpl"): > # Set path to example subdirectory: > path:="examples\/": > # rand_path:="debug_exs\/": # debugging the spreadsheet > rand_path:=path: # random examples in paper > > # Factorization improvement routines: > read("improve_factorization-3var.mpl"): > > infolevel[imfac]:=1; > Digits:=floor(evalhf(Digits)); > 1 14 > # A note about backward_error(F, G, [vars]): > # We compute c so that ||F-c*G||_2 is as small as possible. The > output is ||F - c*G||_2 / ||F||_2. > > # Example 15a: > read(cat(rand_path,"exF15aK2K")): > F15a := expand(evalf(cleanF15a + noiseF15a));: > norm(evalf(noiseF15a),2) / norm(expand(F15a),2); > ceil(log10(%)); > 4 4 4 2 2 2 81. x + 16. y - 648.00300000000 z + 72. x y - 648. x 2 2 2 - 288. y + 1296. + 0.0020000000000000 x z 2 2 2 + 0.0010000000000000 y z - 0.0070000000000000 z -5 0.49089262694426 10 -5 > # First Round (SVD) of factorization > #read(cat(rand_path,"exF15a-factors")); > > t15a:=time(): > facF15a:=appfac(F15a,[x,y,z]); > print('number_of_factors'=nops(facF15a)); > t15a:=time()-t15a; > c,berror15a:=backward_error(F15a, convert(facF15a, `*`), [x,y,z]); > > ceil(log10(berror15a)); > ` the biggest gap, the last r-th singular values and the number of factors `, 1523328000660.7, 0.17194824923987e-13, 2 ` The time for computing the number of factors*****`, 2.734 `The time for the entire factorization******`, 6.319 2 -15 [0.748440129241907614 z + 0.136232381504041622 10 y z 2 -15 + 0.117605170424153624 y + 0.505928041775944784 10 x z -15 2 + 0.162861561433393287 10 x y + 0.264611777828068606 x -15 - 0.415686352202264540 10 z -15 - 0.339914217312790158 10 y -17 + 0.366118562152615828 10 x - 1.05844768880722740, 2 -15 0.748435931568262936 z + 0.136232806701549390 10 y z 2 -15 - 0.117605665816576346 y - 0.625057607547666410 10 x z -15 2 + 0.228417822717985706 10 x y - 0.264612603713145656 x -15 - 0.415687649608930032 10 z -15 - 0.339915278225587526 10 y -17 - 0.976300165280065290 10 x + 1.05844983735575670] number_of_factors = 2 6.359 -12 -1156.8187577721, 0.50423373373502 10 -12 > evalf(facF15a,10); 2 -15 2 [0.7484401292 z + 0.1362323815 10 y z + 0.1176051704 y -15 -15 + 0.5059280418 10 x z + 0.1628615614 10 x y 2 -15 + 0.2646117778 x - 0.4156863522 10 z -15 -17 - 0.3399142173 10 y + 0.3661185622 10 x - 1.058447689, 2 -15 2 0.7484359316 z + 0.1362328067 10 y z - 0.1176056658 y -15 -15 - 0.6250576075 10 x z + 0.2284178227 10 x y 2 -15 - 0.2646126037 x - 0.4156876496 10 z -15 -17 - 0.3399152782 10 y - 0.9763001653 10 x + 1.058449837] > Digits:=floor(evalhf(Digits)); 14 > # Improvement Round (Gauss-Newton) of Factorization > > t15an:=time(): > > facF15a:=improve_factorization(F15a,facF15a,[x,y,z]);: > t15an:=time()-t15an + t15a; > c,berror15a:=backward_error(F15a, convert(facF15a, `*`), [x,y,z]); > ceil(log10(berror15a)); > improve_factorization: -12 org_error = 0.50423373373502 10 improve_factorization: iters = 2 improve_factorization: -13 new_error = 0.31377644409897 10 , improvement = 16.069840270609 -28 -28 [35.999963461762 - 0.404449990787 10 z - 0.242909794654 10 y -28 2 - 0.1426413252199 10 x - 25.455974434027 z -28 2 + 0.192710329307 10 y z - 3.9999915753538 y -28 -28 + 0.188623566552 10 x z - 0.242003648332 10 x y 2 - 8.9999859549935 x , 36.000036538274 -28 -29 - 0.108040873630 10 z + 0.92329072118 10 y -28 2 - 0.2808784934930 10 x + 25.455831662599 z -29 2 + 0.37670865379 10 y z - 4.0000084246596 y -29 -28 - 0.44777907289 10 x z + 0.169533193867 10 x y 2 - 9.0000140450262 x ] 6.880 -13 0.99999999999992, 0.77367637506122 10 -13 > > save facF15a, t15an, berror15a, cat(rand_path,"exF15a-factors-im"); > > facF15a; -28 -28 [35.999963461762 - 0.404449990787 10 z - 0.242909794654 10 y -28 2 - 0.1426413252199 10 x - 25.455974434027 z -28 2 + 0.192710329307 10 y z - 3.9999915753538 y -28 -28 + 0.188623566552 10 x z - 0.242003648332 10 x y 2 - 8.9999859549935 x , 36.000036538274 -28 -29 - 0.108040873630 10 z + 0.92329072118 10 y -28 2 - 0.2808784934930 10 x + 25.455831662599 z -29 2 + 0.37670865379 10 y z - 4.0000084246596 y -29 -28 - 0.44777907289 10 x z + 0.169533193867 10 x y 2 - 9.0000140450262 x ] > cleanF15a; 4 4 4 2 2 2 2 81 x + 16 y - 648 z + 72 x y - 648 x - 288 y + 1296 > alias(sqrt2=RootOf(z^2-2)): facs:=evala(Factor(cleanF15a,sqrt2)); 2 2 2 16 (y + 9/4 x - 9 - 9/2 sqrt2 z ) 2 2 2 (y - 9 + 9/4 x + 9/2 sqrt2 z ) > f1 := evalf(4*op(2,facs)); 2 2 2 4. y + 9. x - 36. - 25.455844122716 z > evalf(facF15a[2]+f1,10); -28 -29 0.00003654 - 0.1080408736 10 z + 0.9232907212 10 y -28 2 - 0.2808784935 10 x - 0.00001246 z -29 -5 2 + 0.3767086538 10 y z - 0.8425 10 y -29 -28 - 0.4477790729 10 x z + 0.1695331939 10 x y 2 - 0.000014045 x # #