\\ Implementation of Method 3 from the Arith 17 \\ article "Correctly rounded multiplication by arbitrary \\ precision constants" \\ \\ print("First choose a precision at least equal to the double of the precision n "); { nearest(z,n,k,a,e)=; e=floor(log(abs(z))/log(2)); a=2^(n-1-e); k=round(a*z); if(abs(z*a-k)-1/2,k/a,if(k%2,(k-1)/a,k/a)); } { ulp(z,n,e,u)=; e=floor(log(abs(z))/log(2)); u=2^(e-n+1) } { expo(z,e)=; e=floor(log(abs(z))/log(2)) } { ch(C,n)=; nearest(C,n) } { cl(C,n,a)=; \\a=z-ch(z,n); nearest(C-ch(C,n),n); } { u1(C,x,n)=; nearest(cl(C,n)*x,n); } { u2(C,x,n)=; nearest(ch(C,n)*x+u1(C,x,n),n); } { echec(C,x,n)=; u2(C,x,n)-nearest(C*x,n) } { method3(const,n,C,Ch,Cl,a,b,Xcut,A,B)=; C=const; /* \p5*n; */ while(C>= 2, C =C/2); while(C<1, C =C*2); Ch=ch(C,n);Cl=cl(C,n);print("Exponent of ulp(Cl) = ",expo(ulp(Cl,n))); A=2*(Ch+Cl);B=2^n*(Ch+Cl); Xcut=floor((2^n)/C); print(); print("2^(n-1) = ",2^(n-1)); print("Xcut = ",Xcut); print("2^(n) = ",2^(n)); print(); print("For x between 1 and xcut"); if(ulp(Cl,n)>=1/2^(2*n+1), print(); print("We solve the first inequality and check if the solutions are bad cases for the algorithm"); solve_ineq(A,2^(n)*ulp(Cl,n), 2^(n+1)*ulp(Cl,n) ,2^(n-1),Xcut,n,C); print("Done"); print(); print("We solve the second inequality and check if the solutions are bad cases for the algorithm"); solve_ineq(A,0,2^(n+1)*ulp(Cl,n),2^(n-1),Xcut,n,C); print("Done"); print(); print("We solve the third inequality and check if the solutions are bad cases for the algorithm"); solve_ineq(A,2^(n+1)*ulp(Cl,n),2^(n+1)*ulp(Cl,n) ,2^(n-1),Xcut,n,C); print("Done"), print("It suffices to check the convergents")); print(); print("For x between xcut and 2"); if(ulp(Cl,n)>=1/2^(2*n), print(); print("We solve the inequality and check if the solutions are bad cases for the algorithm"); solve_ineq(A/2,2^(-n-1),2^(-n),Xcut+1,2^n-1,n,C); print("Done"); , print("It suffices to check the convergents")); } { solve_ineq(a,b,borne,Xdepart,Xarrivee,n,C,Xdep,Xarr,Liste)=; Xarr=Xarrivee; Xdep=Xdepart; while(Xarr>=Xdep, Liste=dist_contrainte(a,b,borne,Xdep,Xarr,n,C); if(Liste[1]>borne, break,); if(Liste[2]>Xarr, break,); if(Liste[3]!=0, \\print(); print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"); print("The algorithm fails for X = ",Liste[2]); print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"); \\print() , \\print("Algo OK") ); Xdep=Liste[2]+1; ); } { dist_contrainte(A,B,borne,Xdepart,Xarrivee,n,C,N,a,b,x,y,d,u,v,U,V,coeff1)=; print(); a=A;b=B+A*Xdepart; x=frac(-a);y=1-frac(-a);d=frac(b);u=1;v=1;coeff1=0; N=Xarrivee-Xdepart; while(u+v<2*N, if(d