> restart: > read("appbifac0412.mpl"): > read("backwards_error.mpl"): > infolevel[appfac]:=1: # > ## Nagasaka ISSAC'02 > f1:=(x^2+u*x+2*u-1); > f2:=(x^3+u^2*x-u+7); > fnoise:=0.2*x; 2 f1 := x + u x + 2 u - 1 3 2 f2 := x + u x - u + 7 fnoise := 0.2 x > fnag:=expand(f1*f2+fnoise); 5 3 2 2 2 4 3 2 2 fnag := x + x u - x u + 7 x + u x + u x - 2 u x + 7 u x 3 3 2 3 + 2 u x + 2 u x - 2 u + 15 u - x - 7 + 0.2 x > > facsnag:=appfac(fnag,[x,u]); > appfac: sv = 0.00168900498079447019 appfac: e0 = 2 appfac: maxrate = 10.66697691 appfac: nd = 2 ` the time for computing the number of factors*****`, .60e-1 `the time for the entire factorization******`, 2.230 appfac: the factorization is correct facsnag := [1.430595079 + 0. I - (0.006771717432 + 0. I) x 2 - (0.2046885565 + 0. I) u - (0.009135194955 + 0. I) x + (0.00825547636704574046 + 0. I) x u 2 3 - (0.006356747617 + 0. I) u + (0.2067944755 + 0. I) x 2 - (0.000380983847633870043 + 0. I) x u 2 + (0.206487613385874208 + 0. I) x u 3 - (0.0003849527731 + 0. I) u , 0.4241319523 + 0. I - (0.01161405484 + 0. I) x - (0.8479885642 + 0. I) u 2 - (0.4231138473 + 0. I) x - (0.427405435436850611 + 0. I) x u 2 - (0.0005657933935 + 0. I) u ] > # Size of the Noise: > > backward_error(f1*f2, f1*f2+fnoise,[x,u]); 0.01007534301 > # Backward error of the computed factorization: > > backward_error(fnag, convert(facsnag,`*`),[x,u]); 0.009454097130 > > # Backwards Error of the Factors: > > backward_error(f1, facsnag[2], [x,u]); 0.01078377962 > backward_error(f2, facsnag[1], [x,u]); 0.01068382715 # # > ## Bigger Random Example: # > myrandpoly := proc(df, range) > randpoly([x,y],degree=df,coeffs=rand(-range..range),dense); > end proc: # > df := 4: fa:= myrandpoly(df, 5); 4 3 3 2 2 2 3 2 fa := 5 - 4 x - 2 x y - 2 x - 4 x y - 4 x + 4 x y + 4 x y 4 3 2 - x y - 5 x - 2 y - 5 y + y + 4 y > dg := 4: ga:= myrandpoly(dg, 5); 3 3 2 2 3 2 ga := -1 - x y - 4 x + x y + 3 x - 3 x y + 4 x y - x y - x 4 3 2 + y - 2 y + 3 y > F:=expand(fa*ga); 2 3 5 4 7 6 8 F := -5 - 4 y + 24 x - 29 x + 14 x + 14 x + 16 x - 4 x - 2 y 7 6 3 2 2 6 2 5 - y + 5 y - 8 x y - 15 x y - 16 x y + 10 x y 2 4 2 3 7 6 5 4 - 12 x y - 3 x y + 10 x y + 3 x y - 18 x y + 8 x y 3 2 2 2 3 2 2 - 23 x y - 24 x y + 23 x y + 8 x y - 4 x y + 14 y 3 4 7 6 5 3 5 2 + 7 y + 2 y + 4 x y + 6 x y + 16 x y - 2 x y 4 4 4 3 4 2 6 2 4 - 2 x y - 18 x y - 45 x y + 2 x y + 9 x y 3 5 3 4 3 3 5 + 12 x y + 3 x y + 49 x y - 13 y > dF := 4: > ua:=1./10^(10+1.5)*myrandpoly(dF,10^10); 4 3 ua := -0.01387302120 + 0.02631739991 x - 0.0004887373488 x y 3 2 2 - 0.01101532645 x + 0.02422310548 x y 2 2 - 0.006748276227 x y + 0.03080729005 x 3 2 + 0.02209525436 x y + 0.02180093499 x y 4 + 0.01107778191 x y - 0.007940033503 x + 0.02046402043 y 3 2 + 0.02531580955 y - 0.002639374096 y - 0.02055921793 y > f:=evalf(expand(F+ua)); 2 f := -0.007940033503 x - 4.020559218 y + 24.03080729 x 3 5 4 7 6 - 29.01101533 x + 14. x + 14.02631740 x + 16. x - 4. x 8 7 6 3 2 - 2. y - 1. y + 5. y - 7.988922218 x y - 15. x y 2 6 2 5 2 4 2 3 7 - 16. x y + 10. x y - 12. x y - 3. x y + 10. x y 6 5 4 3 + 3. x y - 18. x y + 8. x y - 23.00048874 x y 2 2 2 3 - 23.97577689 x y + 22.99325172 x y + 8.022095254 x y 2 2 - 3.978199065 x y - 5.013873021 + 13.99736063 y 3 4 7 6 + 7.025315810 y + 2.020464020 y + 4. x y + 6. x y 5 3 5 2 4 4 4 3 4 2 + 16. x y - 2. x y - 2. x y - 18. x y - 45. x y 6 2 4 3 5 3 4 3 3 + 2. x y + 9. x y + 12. x y + 3. x y + 49. x y 5 - 13. y # > fac:=appfac(f,[x,y]); appfac: sv = 0.000351199177146201743 appfac: e0 = 3 appfac: maxrate = 239.8343791 appfac: nd = 2 ` the time for computing the number of factors*****`, .370 `the time for the entire factorization******`, 5.360 appfac: the factorization is correct fac := [0.5645723523 + 0. I + (0.4519894562 + 0. I) y 2 - (0.5637167618 + 0. I) x + (0.1135659113 + 0. I) y - (0.111690940810708364 + 0. I) y x 2 3 - (0.4507679843 + 0. I) x - (0.5643339274 + 0. I) y 2 + (0.451314301018682973 + 0. I) y x 2 + (0.000212967895177447681 + 0. I) y x 3 4 - (0.2257125206 + 0. I) x - (0.2257899663 + 0. I) y 3 + (0.451707451823876283 + 0. I) y x 2 2 - (0.450837317536006421 + 0. I) y x 3 - (0.225227377559899439 + 0. I) y x 4 - (0.4513263856 + 0. I) x , -0.1815304323 + 0. I - (0.0002753749700 + 0. I) y - (0.1819834964 + 0. I) x 2 + (0.5438561420 + 0. I) y - (0.181344423097735352 + 0. I) y x 2 3 + (0.5439093800 + 0. I) x - (0.3624291521 + 0. I) y 2 + (0.725333891306389477 + 0. I) y x 2 + (0.180628395317499357 + 0. I) y x 3 4 - (0.7258674964 + 0. I) x + (0.1810907592 + 0. I) y 3 - (0.544294630135757518 + 0. I) y x 2 2 + (0.0000864379990762651715 + 0. I) y x 3 - (0.181386556824595452 + 0. I) y x 4 - (0.0001163437293 + 0. I) x ] > nops(fac); 2 > ftilde := expand(convert(fac,`*`)): > backward_error(F,f,[x,y]); > round(log10(%)); > 0.0006863719221 -3 > backward_error(f,ftilde,[x,y]); > round(log10(%)); 0.0005755697678 -3 > # # # # # # #