> restart;
> read "multifac_1.2.mpl":
> path:="examples\/":
# Zeng's sixth example
> p1:=x^2+x*y-y^2+1;

                              2          2
                       p1 := x  + x y - y  + 1

> p2:=expand((x^3+1)*(y^2-2));

                            3  2      3    2
                     p2 := x  y  - 2 x  + y  - 2

> p3:=x^3-y^3-3*x*y^2+2;

                             3    3        2
                      p3 := x  - y  - 3 x y  + 2

> pp:=randpoly([x,y],degree=4,coeffs=rand(-5..5),terms=10);

                      2              2        3      4      3
  pp := -4 - 2 y + 3 x  + 4 x y + 4 x  y - 3 y  - 3 x  + 4 x  y

              3
         - x y

> cleanF14:=p1^4*p2^3*p3:
> noiseF14:='`*`'(10^(-8),pp):
> save cleanF14, noiseF14, cat(path,"exF14");
> zengex6:=evalf(expand(cleanF14+noiseF14)):
# The example has univariate  multiple factors  (x^3+1)^3, (y^2-2)^3. 
> t14:=time():
> res:=multifac_cont(zengex6,[x,y],8);
> t14:=time()-t14;
` the biggest gap, the last r-th singular values and the number of
factors `, 19325210.41, .6474953137e-9, 5
` The time for computing the number of factors*****`, .211
`The time for the entire factorization******`, 17.986

                                               6
  res := [[5.649026103 - 0.706128264342899548 y

                                  -8  5                        4
         + 0.128501856795345096 10   y  + 4.23676957390788012 y

                                  -8  3                        2
         + 0.253707824588418993 10   y  - 8.47353916157238630 y

                                  -8
         + 0.190331440121101876 10   y, 1], [

                                        3
        -(0.841421144403573407 + 0. I) y

                                           2
         - (2.52426344349612331 + 0. I) x y

                                   -8          2
         + (0.462709134174881798 10   + 0. I) x  y

                                          3
         + (0.841421147695857075 + 0. I) x

                                   -8          2
         - (0.173980778590697766 10   + 0. I) y

                                   -8
         + (0.640162425863826717 10   + 0. I) x y

                                   -9          2
         + (0.677962755348901756 10   + 0. I) x

                                   -10
         + (0.865744232832211152 10    + 0. I) y

                                   -8
         - (0.594832540312230908 10   + 0. I) x + 1.682842295 + 0. I,

        1], [

                                 -8                          -8
        -(0.124393252924683672 10   - 0.228678581120713224 10   I) y

         + (0.676200707702497850 - 1.34979274479615530 I) x

         + 0.8308544429 + 1.260503361 I, 3], [

                                -10                          -8
        (0.316390240141101752 10    + 0.260250315219718716 10   I) y

         + (0.0332721122190638985 - 1.50933124818121067 I) x

         - 1.323755250 + 0.7258511324 I, 3], [

                                 -8
        -(0.707478435280108376 10   + 0. I) y

         - (3.57884843402732588 + 0. I) x - 3.578848432 + 0. I, 3], [

                                      2
        (1.07457859336006866 + 0. I) y

         - (1.07457837253731126 + 0. I) x y

                                         2
         - (1.07457860806282324 + 0. I) x

                                   -6
         - (0.107178052399655740 10   + 0. I) y

                                   -6
         + (0.201806598217219931 10   + 0. I) x - 1.074578529 + 0. I,

        4]]


                            t14 := 73.516

> result:=expand(res[1][1]*res[2][1]*res[3][1]^3*res[4][1]^3*res[5][1]^3
> *res[6][1]^4):
> backward_error(result,zengex6,[x,y]);
Approximate factorization is not a real polynomial

                   c = -47.80324096 - 427.2862033 I


                                         -6
                          0.2087974882 10

> h1:=res[6][1];

                                      2
  h1 := (1.07457859336006866 + 0. I) y

         - (1.07457837253731126 + 0. I) x y

                                         2
         - (1.07457860806282324 + 0. I) x

                                   -6
         - (0.107178052399655740 10   + 0. I) y

                                   -6
         + (0.201806598217219931 10   + 0. I) x - 1.074578529 + 0. I

> backward_error(h1,p1,[x,y]);
Approximate factorization is not a real polynomial

                           c = -1.074578526


                                         -6
                          0.1371203304 10

> h2:=expand(res[1][1]*res[3][1]^3*res[4][1]^3*res[5][1]^3):
> backward_error(h2,expand(p2^3),[x,y]);
Approximate factorization is not a real polynomial

                   c = -42.60801611 - 380.8490289 I


                                         -7
                          0.2370788088 10

> h3:=res[2][1];

                                        3
  h3 := -(0.841421144403573407 + 0. I) y

                                           2
         - (2.52426344349612331 + 0. I) x y

                                   -8          2
         + (0.462709134174881798 10   + 0. I) x  y

                                          3
         + (0.841421147695857075 + 0. I) x

                                   -8          2
         - (0.173980778590697766 10   + 0. I) y

                                   -8
         + (0.640162425863826717 10   + 0. I) x y

                                   -9          2
         + (0.677962755348901756 10   + 0. I) x

                                   -10
         + (0.865744232832211152 10    + 0. I) y

                                   -8
         - (0.594832540312230908 10   + 0. I) x + 1.682842295 + 0. I

> backward_error(h3,p3,[x,y]);
Approximate factorization is not a real polynomial

                           c = 0.8414211473


                                         -8
                          0.3230606904 10

> facF14:=res:
> save facF14, t14, cat(path,"exF14-factors");
> 
> 
>