getdigit := proc(hatw);if hatw > 0.5 then d := 1 elif hatw < -0.5 then d := -1 else d := 0 end if;d;end;with(RandomTools);oneiteration_polynom := proc(wr,wi,x,y,n,delta,Bigdr,Bigdi);for k from 0 to n do
# we simulate the fact that hatwr is an approximation to hatwr
# at a distance bounded by delta/2
epsilon := Generate(float(range=0..delta/2,digits=5)); hatwr[k] := wr[k]+epsilon;
epsilon := Generate(float(range=0..delta/2,digits=5)); hatwi[k] := wi[k]+epsilon;
dr[k] := getdigit(hatwr[k]);
NewBigdr[k] := 2*Bigdr[k]+dr[k];
di[k] := getdigit(hatwi[k]);
NewBigdi[k] := 2*Bigdi[k]+di[k];end do;# now we update wr and wifor k from 0 to n-1 do newwr[k] := 2*(wr[k]-dr[k]+x*dr[k+1]-y*di[k+1]); newwi[k] := 2*(wi[k]-di[k]+y*dr[k+1]+x*di[k+1]);end do;newwr[n] := 2*(wr[n]-dr[n]);newwi[n] := 2*(wi[n]-di[n]);(newwr,newwi,NewBigdr,NewBigdi);end;complexEmethod_polynom := proc(p::polynom,t::complex, delta, niter);n := degree(p,z);
x := Re(t);
y := Im(t);for k from 0 to n do wr[k] := Re(coeff(p,z,k));
wi[k] := Im(coeff(p,z,k));
Bigdr[k] := 0;
Bigdi[k] := 0;
od;for j from 1 to niter do
(wr,wi,Bigdr,Bigdi) := oneiteration_polynom(wr,wi,x,y,n,delta,Bigdr,Bigdi);
# printf("ITERATION %a\n",j);
# for k from 0 to n do printf("wr[%a] = %a\n",k,wr[k]) end do;
# for k from 0 to n do printf("wi[%a] = %a\n",k,wi[k]) end do;end do;realpart := Bigdr[0]/2^(niter-1);imaginarypart := Bigdi[0]/2^(niter-1);
sol := realpart +I*imaginarypart;sol;end;p := (1+I)*z^3-(0.5+1.25*I)*z^2+z+1;funcp :=z -> (1+I)*z^3-(0.5+1.25*I)*z^2+z+1;complexEmethod_polynom(p,0.01+I*0.1,0.00000001,20);evalf(%);evalf(funcp(0.01+I*0.1));Digits := 45;evalf(%);p;latex(1.018121+.110106*I);oneiteration_rational := proc(wr,wi,qr,qi,x,y,n,delta,Bigdr,Bigdi);for k from 0 to n do
# we simulate the fact that hatwr is an approximation to hatwr
# at a distance bounded by delta/2
epsilon := Generate(float(range=0..delta/2,digits=5)); hatwr[k] := wr[k]+epsilon;
epsilon := Generate(float(range=0..delta/2,digits=5)); hatwi[k] := wi[k]+epsilon;
dr[k] := getdigit(hatwr[k]);
NewBigdr[k] := 2*Bigdr[k]+dr[k];
di[k] := getdigit(hatwi[k]);
NewBigdi[k] := 2*Bigdi[k]+di[k];end do;# now we update wr and wi
newwr[0] := 2*(wr[0]-dr[0]+x*dr[1]-y*di[1]);
newwi[0] := 2*(wi[0]-di[0]+y*dr[1]+x*di[1]);for k from 1 to n-1 do newwr[k] := 2*(wr[k]-dr[k]-qr[k]*dr[0] + qi[k]*di[0]+x*dr[k+1]-y*di[k+1]); newwi[k] := 2*(wi[k]-di[k]-qi[k]*dr[0] - qr[k]*di[0]+y*dr[k+1]+x*di[k+1]);end do;newwr[n] := 2*(wr[n]-dr[n]-qr[n]*dr[0] + qi[n]*di[0]);newwi[n] := 2*(wi[n]-di[n]-qi[n]*dr[0] - qr[n]*di[0]);(newwr,newwi,NewBigdr,NewBigdi);end;complexEmethod_rational := proc(numerator::polynom,denominator::polynom,t::complex, delta, niter);n := degree(numerator,z);
x := Re(t);
y := Im(t);for k from 0 to n do wr[k] := Re(coeff(numerator,z,k));
wi[k] := Im(coeff(numerator,z,k));
qr[k] := Re(coeff(denominator,z,k));
qi[k] := Im(coeff(denominator,z,k));
Bigdr[k] := 0;
Bigdi[k] := 0;
od;for j from 1 to niter do
(wr,wi,Bigdr,Bigdi) := oneiteration_rational(wr,wi,qr,qi,x,y,n,delta,Bigdr,Bigdi);
# printf("ITERATION %a\n",j);
# for k from 0 to n do printf("wr[%a] = %a\n",k,wr[k]) end do;
# for k from 0 to n do printf("wi[%a] = %a\n",k,wi[k]) end do;end do;realpart := Bigdr[0]/2^(niter-1);imaginarypart := Bigdi[0]/2^(niter-1);
sol := realpart +I*imaginarypart;sol;end;pnum := (1+1.3*I)*z^4 + 0.7*z^3 + I*z^2 + (0.5-I)*z + (1-I);pden := 1 + (0.06 - 0.03*I)*z + 0.08*z^2 - (0.04+0.03*I)*z^3 + 0.07*I*z^4;f := z -> ((1+1.3*I)*z^4 + 0.7*z^3 + I*z^2 + (0.5-I)*z + (1-I))/(1 + (0.06 - 0.03*I)*z + 0.08*z^2 - (0.04+0.03*I)*z^3 + 0.07*I*z^4);t := 0.05 - 0.02*I;complexEmethod_rational(pnum,pden,t, 0.125, 30);evalf(%);f(t);