"Lemma 2.2: bounds for the function g given by    g(a,b,c,x) = kappa(a/(a^(2)\

  +x^(2))+(b)/(b^(2)+x^(2)))-(c)/(c^(2)+x^(2))  kappa=1/(sqrt(5)),  

   a[0]=(sqrt(5)-1)/(2) , b[0]=(sqrt(5)+1)/(2) , c[0]=1 "
> restart;
> 
> #part (i)
> kappa := 1/sqrt(5):
> a0 := (sqrt(5)-1)/2:
> b0 := (sqrt(5)+1)/2:
> c0 := 1:
> x0 := fsolve( 2*x^6+4*x^4-1 , x=0..infinity); # root of the polynomial 2x^6+4x^4-1
> g := proc(a,b,c,x) kappa*(a/(a^2+x^2)+b/(b^2+x^2))-c/(c^2+x^2) end:
> simplify(g(a0,b0,c0,0)) ;
> g0:= evalf(g(a0,b0,c0,x0)) ;
                                0.6720163413
                                      0
                                -0.1215851071
"Application of Lemma 2.2 to get Lemma 2.3.   epsilon=0.15, 1<sigma<1+epsilon\

  , 1-`&epsilon;`<beta<1,    sigma[1](x) =(1+sqrt(1+4 x^(2)) )/(2)    

   a=sigma[1]-beta, b=sigma[1]-1+ beta, c=sigma-1+ beta, and x=|k-1| 

   gamma[0],   m[1],m[2],m[3] are given by |a-a[0]| <m[1 ]`&epsilon;`,  

   |b-b[0]| <m[2 ]`&epsilon;`,  |c-c[0]| <m[3 ]`&epsilon;`,  (a-a[0] 

   )/(`&epsilon;`) <(sigma[1](1+`&epsilon;`)-(1-`&epsilon;`)-a[0]  

   )/(`&epsilon;`) , (b-b[0] )/(`&epsilon;`) <(sigma[1](1+`&epsilon;`)-1+1  

   )/(`&epsilon;`) , (c-c[0] )/(`&epsilon;`) <((1+`&epsilon;`)-1+1-1  

   )/(`&epsilon;`)   m[0] = (kappa m[1])/((a[0]-m[1]eps)^2)+ (kappa 

   m[2])/((b[0]-m[2]eps)^2)+(m[3])/((c[0]-m[3]eps)^2)"
> 

> #p9-10 - #proof of Lemma 2.3.
> eps:=0.15:
> s1:=proc(x) ((1+sqrt(1+4*x^2))*(1/2)) end:
> 
> m1 := evalf( ( s1(1+eps) - (1-eps) - a0 ) / eps) ;
> m2 := evalf( ( s1(1+eps) - b0 ) / eps) ;
> m3 := 1;
> m0:= evalf( kappa*m1/(a0-m1*eps)^2+kappa*m2/(b0-m2*eps)^2+m3/(c0-m3*eps)^2);
> alpha1 := evalf(m0*eps);
                                 1.906397542
                                 0.906397542
                                      1
                                 9.300018076
                                 1.395002711
> #p10-11 - #proof of Lemma 2.4.
> restart:
> eps:=0.15:
> kappa := 1/sqrt(5):
> s1:=proc(x) ((1+sqrt(1+4*x^2))*(1/2)) end:
> 
> #finding bound alpha20 for 1/s - kappa/s1 - kappa/(s1-1)
> f :=proc(x) 1/x - kappa/s1(x) - kappa/(s1(x)-1) end:
> maximize( f(x) ,x=1..1+eps , location);
> alpha20 := evalf(f(1+eps));
> 
> g := proc(a,b,c,x) kappa*(a/(a^2+x^2)+b/(b^2+x^2))-c/(c^2+x^2) end:
> a0 := (sqrt(5)-1)/2:
> b0 := (sqrt(5)+1)/2:
> c0 := 1:
> x0 := fsolve( 2*x^6+4*x^4-1 , x=0..infinity): 
> g0:= evalf(g(a0,b0,c0,x0)) :
> m1 := evalf( ( s1(1+eps) - (1-eps) - a0 ) / eps) :
> m2 := evalf( ( s1(1+eps) - b0 ) / eps) :
> m3 := 1:
> m0:= evalf( kappa*m1/(a0-m1*eps)^2+kappa*m2/(b0-m2*eps)^2+m3/(c0-m3*eps)^2):
> 
> #finding bound alpha21 for g(sigma1-1,sigma1,sigma,kgamma0)
> alpha21:=-evalf(g0-m0*eps);
> 
> #finding bound alpha22 for g(sigma1-1,sigma1,sigma,kgamma0)
> alpha22:=evalf(alpha21+eps);
> 
> 
                 0.0214699497, {[{x = 1.15}, 0.0214699497]}
                                0.0214699497
                                 1.516587818
                                 1.666587818
> 
> restart;
> kappa:=evalf(1/sqrt(5)): 
> evalf(1-kappa)/2:
> print("(1-kappa)/2"=%);
                        "(1-kappa)/2" = 0.2763932023
                           Case1; >(gamma[0], 1);
> restart;
> kappa:=evalf(1/sqrt(5)): 
> 
> # values for alpha_i's (Lemma 2.3-2.4)
> alpha20 := 0.021470:
> alpha21 := 1.516588:
> alpha22 := 1.666588:
> 
> # values for d(0) and D(k) (Lemma 2.5)
> d0 := -0.0512 :
> vD[1] := 0.3918 :
> vD[2] := 0.3915 :
> vD[3] := 0.4062 :
> vD[4] := 0.4266 :
> 
> # values for trigonometric polynomial
> va:=0.8924:
> vb:=0.1768:
> Q:=8*(va+cos(t))^2*(vb+cos(t))^2: # definition of the trigonometric polynomial with respect to va, b
> P:=combine(Q,trig): # linear form of the polynomial to obtain a0,a1,a2,a3,a4  
> a[1]:=evalf(coeff(P,cos(t))):
> a[2]:=evalf(coeff(P,cos(2*t))):
> a[3]:=evalf(coeff(P,cos(3*t))):
> a[4]:=evalf(coeff(P,cos(4*t))):
> a[0]:=combine( P-(a[1]*cos(t)+a[2]*cos(2*t)+a[3]*cos(3*t)+a[4]*cos(4*t)) ,trig):
> print(["a0"= a[0],"a1"= a[1],"a2"= a[2],"a3"= a[3],"a4"= a[4]]);
> 
> 
> # value for r as given p13
> r:=evalf(sqrt(a[0])/(sqrt(a[1])-sqrt(a[0]))):
> print("r"= r);
> 
>  #coefficient of (log d_K)
> c1:=evalf( (sum(a[k], k = 0 .. 4) ) * ((1-kappa)*(1/2)) ) :
> print("c_1"= c1);
> evalf( c1/( a[1]/(1+r)-a[0]/r ) ) :
> print("coefficient of (log d_K)"= %);
> 
> #coefficient of (log gamma_0)n_K
> c2:=evalf( (sum(a[k], k = 1 .. 4) ) * ((1-kappa)*(1/2)) ) : 
> print("c_2"= c2);
> evalf( c2/( a[1]/(1+r)-a[0]/r ) ) :
> print("coefficient of (log gamma_0)n_K"= %);
> 
> #coefficient of nK
> c3:=evalf( a[0]*(-(1-kappa)/2*log(Pi) +d0) + sum(a[k] *((1-kappa)/2*log(k/Pi)+vD[k]) , k = 1 .. 4) ) :
> print("c_3"= c3);
> evalf( c3/( a[1]/(1+r)-a[0]/r ) ) : 
> print("coefficient of n_K"= %);
> 
> #constant term
> c4:=evalf( a[0]*alpha20 + sum(a[k] , k = 1 .. 4)*alpha22 ) : 
> print("c_4"= c4);
> evalf( c4/( a[1]/(1+r)-a[0]/r ) ) :
> print("constant term"= %);
        ["a0" = 9.034112058, "a1" = 15.52951106, "a2" = 9.834965120, 

          "a3" = 4.276800000, "a4" = 1.]
                              "r" = 3.214389997
                             "c_1" = 10.96600761
                  "coefficient of (log d_K)" = 12.54180965
                             "c_2" = 8.469040446
               "coefficient of (log gamma_0)n_K" = 9.686031323
                             "c_3" = 2.649026624
                     "coefficient of n_K" = 3.029688548
                             "c_4" = 51.26034558
                        "constant term" = 58.62639529
> 
                      "Case 2:  (d[2])/(L)<gamma[0] <1"
> restart;
> kappa:=evalf(1/sqrt(5)): 
> 
> # values for d(k) (Lemma 2.5)
> d0 := -0.0512 :
> d[1] := -0.0390 :
> d[2] := 0.2469 :
> d[3] := 0.4452 :
> d[4] := 0.5842 :
> 
> # values for trigonometric polynomial in Case 2
> va:=0.8918:
> vb:=0.1732:
> Q:=8*(va+cos(t))^2*(vb+cos(t))^2: # definition of the trigonometric polynomial with respect to va, b
> P:=combine(Q,trig): # linear form of the polynomial to obtain a0,a1,a2,a3,a4  
> a[1]:=evalf(coeff(P,cos(t))):
> a[2]:=evalf(coeff(P,cos(2*t))):
> a[3]:=evalf(coeff(P,cos(3*t))):
> a[4]:=evalf(coeff(P,cos(4*t))):
> a[0]:=combine( P-(a[1]*cos(t)+a[2]*cos(2*t)+a[3]*cos(3*t)+a[4]*cos(4*t)) ,trig):
> print(["a0"= a[0],"a1"= a[1],"a2"= a[2],"a3"= a[3],"a4"= a[4]]);
> 
> #coeeficient of n_K is negative
> evalf(a[0] *( -(1-kappa)/2*log(Pi) + d0 ) + sum(a[k] *( -(1-kappa)/2*log(Pi) + d[k] ) ,k=1..4)):
> print( % < 0);
> 
> # table values in Case 2
> d2:=2.374:
> r2:=0.248:
> c2 := 1/fsolve( a[0]/r2 - a[1]/(r2+c) + a[1]*r2/(r2^2+d2^2) - a[0]*(r2+c)/((r2+c)^2+d2^2) + sum(a[k],k=0..4)*(1-kappa)/2 ,c=0.0001..1):
> print(["d2"= d2,"r"= r2,"1/c"= c2]);
> 
> # verification of (3.2) and proof of eq (3.4)
> print(d2 > sqrt(r2*(r2+c2))/2); # see eq(3.2) 
> 
> 
> # verification of (3.1) and proof of eq (3.5)
> print(a[0]/(a[1]-a[0])*1/c2 < r2);
> q:= evalf(a[1]/a[0]);
> print(2*a[1] > a[0]);
        ["a0" = 8.963440620, "a1" = 15.41199431, "a2" = 9.772578080, 

          "a3" = 4.260000000, "a4" = 1.]
                              -8.634914714 < 0
              ["d2" = 2.374, "r" = 0.248, "1/c" = 12.73044822]
                            0.8970305400 < 2.374
                            0.1091864462 < 0.248
                                 1.719428394
                          8.963440620 < 30.82398862


                 "Case 3:  (d[1])/(L)<gamma[0] <(d[2])/(L)"
> restart;
> kappa:=evalf(1/sqrt(5)): 
> 
> # values for d(k) (Lemma 2.5)
> d0 := -0.0512 :
> d[1] := -0.0390 :
> d[2] := 0.2469 :
> d[3] := 0.4452 :
> d[4] := 0.5842 :
> 
> # values for trigonometric polynomial in Case 3
> va:=0.8924:
> vb:=0.1771:
> Q:=8*(va+cos(t))^2*(vb+cos(t))^2: # definition of the trigonometric polynomial with respect to va, b
> P:=combine(Q,trig): # linear form of the polynomial to obtain a0,a1,a2,a3,a4  
> a[1]:=evalf(coeff(P,cos(t))):
> a[2]:=evalf(coeff(P,cos(2*t))):
> a[3]:=evalf(coeff(P,cos(3*t))):
> a[4]:=evalf(coeff(P,cos(4*t))):
> a[0]:=combine( P-(a[1]*cos(t)+a[2]*cos(2*t)+a[3]*cos(3*t)+a[4]*cos(4*t)) ,trig):
> print(["a0"= a[0],"a1"= a[1],"a2"= a[2],"a3"= a[3],"a4"= a[4]]);
> 
> # table values in Case 3
> d1:=1.021:
> d2:=2.374:
> r12:=0.236:
> c12 := 1/fsolve( a[0]/r12 - a[1]/(r12+c) + a[1]*r12/(r12^2+d1^2) - a[0]*(r12+c)/((r12+c)^2+d1^2) - a[0]*(r12+c)/((r12+c)^2+d2^2) - a[1]*(r12+c)/((r12+c)^2+4*d2^2) + sum( a[k]*( r12/(r12^2+(k*d1)^2) - (r12+c)/((r12+c)^2+((k-1)*d2)^2) - (r12+c)/((r12+c)^2+((k+1)*d2)^2) ) ,k=2..4) + sum( a[k] ,k=0..4)*(1-kappa)/2 ,c=0.0001..1):
> print(["d1"= d1,"d2"= d2,"r"= r12,"1/c"= c12]);
> 
> #verification of eq (3.7)
> print(a[0]/(a[1]-a[0])*1/c12 < r12);
> 
> #coeeficient of n_K is negative
> evalf( a[0] *( -(1-kappa)/2*log(Pi) + d0 ) + sum(a[k] *( -(1-kappa)/2*log(Pi) + d[k] ) ,k=1..4) ):
> print( %<0);
        ["a0" = 9.039496669, "a1" = 15.53844962, "a2" = 9.839673320, 

          "a3" = 4.278000000, "a4" = 1.]
       ["d1" = 1.021, "d2" = 2.374, "r" = 0.236, "1/c" = 12.73005948]
                            0.1092623214 < 0.236
                              -8.710158576 < 0
                      "Case 4:  0<gamma[0] <(d[1])/(L)"

> restart;
> kappa:=evalf(1/sqrt(5)): 
> 
> # values for d(k) (Lemma 2.5)
> d0 := -0.0512 :
> d[1] := -0.0390 :
> d[2] := 0.2469 :
> d[3] := 0.4452 :
> d[4] := 0.5842 :
> 
> #coeeficient of n_K is negative
> evalf( -(1-kappa)/2*log(Pi) + d0  );
> 
> # trigonometric polynomial in Case 4 is just 1 + cos(theta)
> 
> # the value of r which optimizes c=c(r) is given by the root of the derivative of c(r) (see p17)
> diff( (-(1-kappa)/2*x^2+sqrt(x^2-d1^2*(1+(1-kappa)/2*x)^2)) / (1+(1-kappa)/2*x) ,x):
> R1:=proc(d1,x) (-(1-kappa)*x+(1/4)*(8*x-4*d1^2+4*d1^2*kappa-2*d1^2*x+4*d1^2*x*kappa-2*d1^2*x*kappa^2)/sqrt(4*x^2-4*d1^2-4*d1^2*x+4*d1^2*x*kappa-d1^2*x^2+2*d1^2*x^2*kappa-d1^2*x^2*kappa^2))/(1+(1/2*(1-kappa))*x)-(-(1/2*(1-kappa))*x^2+(1/2)*sqrt(4*x^2-4*d1^2-4*d1^2*x+4*d1^2*x*kappa-d1^2*x^2+2*d1^2*x^2*kappa-d1^2*x^2*kappa^2))*(1/2-(1/2)*kappa)/(1+(1/2*(1-kappa))*x)^2 end:
> 
> # table values in Case 4
> d1:=1.021:
> r1:=fsolve(R1(d1,x),x=0..10):
> c1 := 1/evalf( (-(1-kappa)/2*r1^2 + sqrt(r1^2-d1^2*(1+(1-kappa)/2*r1)^2) ) / (1+(1-kappa)/2*r1) ):
> print(["d1"= d1,"r"= r1,"1/c"= c1]);
> 
> d1:=0.25:
> r1:=fsolve(R1(d1,x),x=0..10):
> c1 := 1/evalf( (-(1-kappa)/2*r1^2 + sqrt(r1^2-d1^2*(1+(1-kappa)/2*r1)^2) ) / (1+(1-kappa)/2*r1) ):
> print(["d1"= d1,"r"= r1,"1/c"= c1]);
> 
> restart;
> kappa:=evalf(1/sqrt(5)): 
> 
> # values for d(k) (Lemma 2.5)
> d0 := -0.0512 :
> d[1] := -0.0390 :
> d[2] := 0.2469 :
> d[3] := 0.4452 :
> d[4] := 0.5842 :
> 
> r1:=1.644:
> d1 := fsolve( (-(1-kappa)/2*r1^2 + sqrt(r1^2-x^2*(1+(1-kappa)/2*r1)^2) ) / (1+(1-kappa)/2*r1)  - x , x=0.01..2):
> c1 := d1:
> print(["1/d1"= 1/d1,"r"= r1,"1/c"= 1/c1]);
                                -0.3675955590
           ["d1" = 1.021, "r" = 2.142611919, "1/c" = 12.54939230]
            ["d1" = 0.25, "r" = 1.534402885, "1/c" = 1.691751135]
          ["1/d1" = 1.999613465, "r" = 1.644, "1/c" = 1.999613465]
                            "Case 5: gamma[0] =0"
> restart;
> kappa:=evalf(1/sqrt(5)): 
> r:=evalf(2*(sqrt(2)-1)/(1-kappa)):
> c1:=evalf( r*(1-(1-kappa)/2*r) / ((1+(1-kappa)/2*r) ) ):
> print(["r"= r,"1/c1"= 1/c1]);
                  ["r" = 1.498638746, "1/c1" = 1.610937637]
>