\\==================================================
\\
\\ PARI/gp commands used to compute and check the computations
\\ in "Explicit versions of smoothed prime ideals theorems ... and applications"
\\ There are four blocks of commands, one for each bound. They can be used either
\\ together or separately.
\\
\\ 2/10/2013 LG & GM
\\
\\==================================================

\\==================================================
\\Upper bound for \sum_{\rho} 1/|\rho(\rho+1)| \leq 0.5375\log\disc_\K - 1.0355n_\K + 5.3879
\\==================================================
\p 231
F(s,g)= 4*(2*s-1)/((2*s-1)^2+4*g^2);
G1(g)=4/sqrt((1+4*g^2)*(9+4*g^2));

G1sq(g)=16/((1+4*g^2)*(9+4*g^2));

s(j)=1+j/2;
H1(q,v,r)=
{
    my(w=i->v^i-v+r,M,V,cr);
    M=matrix(2*q,2*q);
    for(j=1,2*q,M[1,j]=F(s(j),0));
    for(i=1,q-1,for(j=1,2*q,M[2*i  ,j]=F(s(j),w(i))));
    for(i=1,q-1,for(j=1,2*q,M[2*i+1,j]=derivnum(x=w(i),F(s(j),x))));
    for(j=1,2*q,M[2*q,j]=j+1);
    V=vectorv(2*q);
    V[1]=G1(0);
    for(i=1,q-1,V[2*i  ]=G1(w(i)));
    for(i=1,q-1,V[2*i+1]=derivnum(x=w(i),G1(x)));
    V[2*q]=1;
    a=matsolve(M,V);
    
    ar1=ceil(a*10^7)/10^7;
    cr = polsturm(numerator(G1sq('x)-(vector(#ar1,j,F(1+j/2,'x))*ar1)^2),0);
    if (cr > 0, warning("F and g are crossing"); return([]););

    ploth(g=0,4*q,sum(j=1,2*q,ar1[j]*F(s(j),g))-G1(g));
    sum(j=1,#ar1,ar1[j]*[1.,2./s(j)+2./(s(j)-1),psi(s(j)/2),psi((s(j)+1)/2)])
}

A1 = H1(40,1.21,1)
\\%95 =
\\[
\\  0.53747129999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999,
\\  5.38789110852054027698941034250603242209935619025981216770026251662068818026605049888758831135238148327,
\\ -0.68386772836824351949896390069227803635140413632521408647051310085040818802455748892991572636899754273,
\\ -0.15679571290396129326185934742710723437639295864086247158843401329279994657859772930747563035534916201
\\]
[A1[1], -A1[1]*log(Pi)+max(A1[3],(A1[3]+A1[4])/2),A1[2]]
\\%223 = [
\\  0.53747129999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999,
\\ -1.03559118053243112219750590904697786008929523206388757927594351700572686826067593249390339614739216370,
\\  5.38789110852054027698941034250603242209935619025981216770026251662068818026605049888758831135238148327
\\ ]
{
ar1=
[
                                       250548071,
                                    -40769390315,
                                   5795175723671,
                                -642251894123528,
                               54218815728127329,
                            -3508878919641771688,
                           177001043449933176447,
                         -7094015596077453633868,
                        230165538494597837675083,
                      -6150059294311314135993327,
                     137429722146979678372903545,
                   -2603418437013270575777900517,
                   42312055609004270243526202076,
                 -596228700498573506460507915379,
                 7352229299660977983271586796428,
               -79991031610893192700264201347849,
               773449451301413812623754497322110,
             -6689469356634480595952166290773419,
             52049469989158830787111938359894141,
           -366215235328303748457085063911062452,
           2340727373433875029268033585654013101,
         -13647569726889979888635481117558851444,
          72856233722227845138595374556506130213,
        -357308755193444577424236430238048816629,
        1614746152052189321203222537039119640756,
       -6742815290893858601169185495146599758906,
       26081651135346560764877059272421175463793,
      -93662343928951334238283477190373026235970,
      312909679670641654646206585298379548969664,
     -974320711195668140488233654168711408974431,
     2832324810202406292456790051806622242980143,
    -7698430960693278611182246394416801692267482,
    19591848109684436395280873748247452691475281,
   -46741161307608105954759866712283645186774375,
   104654143256889695138455737518470291722254806,
  -220128737522779177621064610160993365216868168,
   435354172671749489946963445292948033362127211,
  -810197596515155479177768566395714272810920262,
  1419759138597775056528221649897613967888797952,
 -2344042914614942938251851695053817440107394201,
  3648003867198618158032688666281279907332926401,
 -5353733754976758827081327735207850440805276490,
  7411481592406340412123547436619012373828347148,
 -9680502044407712819603502062322026576945648999,
 11931531864054985817793303920405879163945577143,
-13877909647596266697860364708436232764310634475,
 15232402120827550086363255671423471322396554451,
-15775334247682258723942059247603410917983570659,
 15412228661977544619641915478149603219261905955,
-14200388071264097711591911264486344481572166054,
 12334252439899208072837355489763427886826472535,
-10094621415908831481370162399779502133799211521,
  7779879804978723458319595088819777701055258725,
 -5642269216814651472704347110867628137752200021,
  3847449822637486914738835007382155275278296082,
 -2464415859390293757851604168538551569779024142,
  1481149204040503957548445332963392835676301519,
  -834221598218683553855012914220968482480130787,
   439683909946941169248931270316639282116138102,
  -216506399095273447319941898397799604643721683,
    99419464389596242263022332671287321846845659,
   -42484842271343000074946235138759717939948915,
    16854882803935701426471290624390658493962917,
    -6191128565930702299565499949693343368179853,
     2099040048795249597742846189024188906401966,
     -654534994786055122946588844807706357840291,
      186947926780001735965997901059271943644593,
      -48675054083574200340006259345314988135384,
       11488195908804597573088392774662830983598,
       -2441545378407438626531756675121759469076,
         463522993226953954733282029125741713819,
         -77845294874645427333933115933295893005,
          11425492896861966116655614216587059464,
          -1443037809175186486864654360574474088,
            153673972446397363248771006965862929,
            -13418974865215897151246028213990153,
               922600572073108333758203469960875,
               -46833786494978206937017577663173,
                 1560648037479364896275100707017,
                  -25610063982827093894391815027
]~/10^7
}
vector(#ar1\2,i,my(k=2*i-1);ceil(solve(x=.0001,10^5,ar1[k]/x^s(k)+ar1[k+1]/x^s(k+1))))
\\%153 = [26479, 12283, 4189, 1607, 714, 359, 199, 119, 75, 50, 34, 25, 18, 13, 10, 8, 6, 5, 4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
S(ar,n)=sum(j=1,#ar,ar[j]/n^s(j));
for(n=2,26500,if(S(ar1,n)<0,print(n)))

\\ for the riemann zeta function:
sum(j=1,#ar1,ar1[j]*(2*derivnum(x=s(j),log(zeta(x)))+log(1/Pi)+2*(1/s(j)+1/(s(j)-1))+psi(s(j)/2)))
\\%7 = 0.0461230560926647732786709783345779368628121452048267272780314366938795204251714504930338499637695
\\ The true value for Riemann zeta function, computed using Odlyzko's list of zeros.
\\ gg=readvec("zeros6.gz");
\\ sum(j=1,max(7*10^5,#gg),G1(gg[j]))*2
\\%20 = 0.0461139618194682210684546716534929680349693203055802951297967834267232767899328107522143
\\ "true" value    0.0461139...
\\ our estimate    0.0461230...


\\==================================================
\\Upper bound for \sum_{\rho} 1/|\rho(\rho+1)(\rho+2)| \leq 0.1763\log\disc_\K - 0.4106 n_\K + 2.2496
\\==================================================
\p 115;
F(s,g)= 4*(2*s-1)/((2*s-1)^2+4*g^2);
G2(g)=8/sqrt((1+4*g^2)*(9+4*g^2)*(25+4*g^2));

G2sq(g)=64/((1+4*g^2)*(9+4*g^2)*(25+4*g^2));

s(j)=1+j/2;
H2(q,r,v)=
{
    my(w=i->v^i-v+r,M,V,cr);
    M=matrix(2*q-1,2*q-1);
    for(j=1,2*q-1,M[1,j]=F(s(j),0));
    for(i=1,q-1,for(j=1,2*q-1,M[2*i  ,j]=F(s(j),w(i))));
    for(i=1,q-1,for(j=1,2*q-1,M[2*i+1,j]=derivnum(x=w(i),F(s(j),x))));
    V=vectorv(2*q-1);
    V[1]=G2(0);
    for(i=1,q-1,V[2*i  ]=G2(w(i)));
    for(i=1,q-1,V[2*i+1]=derivnum(x=w(i),G2(x)));
    a=matsolve(M,V);
    
    ar2=ceil(a*10^7)/10^7;
    cr = polsturm(numerator(G2sq('x)-(vector(#ar2,j,F(1+j/2,'x))*ar2)^2),0);
    if (cr > 0, warning("F and g are crossing"); return([]););
    
    ploth(g=0,4*q,sum(j=1,2*q-1,ar2[j]*F(s(j),g))-G2(g));
    sum(j=1,#ar2,ar2[j]*[1.,2./s(j)+2./(s(j)-1),psi(s(j)/2),psi((s(j)+1)/2)])
}
A2=H2(20,0.75,1.27)
\\%170 =
\\[
\\  0.1762991999999999999999999999999999999999999999999999999999999999999999999999999999999999999993814592388465,
\\  2.249549802190926138340756696020168896514592442098626867510581928264508614007167485626985333716048683882126,
\\ -0.3130595139992637335950866057114697840189112066188633239044252000835179468635905337912078344820712777391851,
\\ -0.10472494174690447166700957436158681271469679231848127686290232004479242672671524469380323703940741705953412
\\]
[A2[1], -A2[1]*log(Pi)+max(A2[3],(A2[3]+A2[4])/2),A2[2]]
\\%223 =
\\[
\\0.1762991999999999999999999999999999999999999999999999999999999999999999999999999999999999999993814592388465,
\\-0.4107071909644246738123950173381914667832527571497913981922582966662358979156318192474944914109839003805493,
\\2.249549802190926138340756696020168896514592442098626867510581928264508614007167485626985333716048683882126
\\]
{
ar2=
[
                      116043280,
                   -15019134746,
                  1306482026256,
                -76315741770330,
               3116274365157230,
             -92621169453588672,
            2074954505670798718,
          -36069656819440696263,
          498302313581120124204,
        -5579712481960141840354,
        51471550420429886034202,
      -396486283111534949768375,
      2579203445202845079404723,
    -14302917461736234191777842,
     68148701393954628073176631,
   -280819396505042268256263139,
   1006203468485334827305158167,
  -3148890161469033145131905085,
   8637410243724442351566255216,
 -20823652528449395097665316823,
  44212581087391037851257051242,
 -82776719893697522350544625956,
 136740298375301487499890195367,
-199275732886794715825307355765,
 255978158207512987528401996503,
-289344568336337774362395707820,
 287063511581435319315875881542,
-249076192247094252611008694964,
 188100860940650555126546470265,
-122861964251233620242612405716,
  68841615651370858094138826161,
 -32737240356857723641726749028,
  13026982181479475895165888915,
  -4255456128181013051676843112,
   1111002720718102002015316745,
   -222834382207437523098078851,
     32228595062589755026085278,
     -2991080884530620994922737,
       133739429590971377317925
]~/10^7
};
vector(#ar2\2,i,my(k=2*i-1);ceil(solve(x=.0001,10^5,ar2[k]/x^s(k)+ar2[k+1]/x^s(k+1))))
\\%153 = [16752, 3413, 884, 303, 126, 60, 31, 17, 10, 6, 4, 3, 2, 1, 1, 1, 1, 1, 1]
S(ar,n)=sum(j=1,#ar,ar[j]/n^s(j));
for(n=1,17000,if(S(ar2,n)<0,print(n)))

\\ for the riemann zeta function:
s(j)=1+j/2;
sum(j=1,#ar2,ar2[j]*(2*derivnum(x=s(j),log(zeta(x)))+log(1/Pi)+2*(1/s(j)+1/(s(j)-1))+psi(s(j)/2)))
\\%6 = 0.001453786587759353798127556395050723616345331560229368055124965207766435084132854589698275111893907072752701
\\ The true value for Riemann zeta function, computed using Odlyzko's list of zeros.
\\ gg=readvec("zeros6.gz");
\\ sum(j=1,max(#gg,max(3*10^5,#gg)),G2(gg[j]))*2
\\%15= 0.0014399631992656256359379806577790068072859743385258053123378097665
\\ the "true" value   0.0014399...
\\ our estimate       0.0014537...



\\==================================================
\\Upper bound for |B_K| = (1/2)\sum_{\gamma} 4/(1+4\gamma^2) \leq 0.5155\log\disc_\K - 1.2432 n_\K + 9.3419
\\==================================================
\p 115
F(s,g)= 4*(2*s-1)/((2*s-1)^2+4*g^2);
G3(g)=2/(1+4*g^2);
s(j)=1+j/2;
H3(q,v)=
{
    my(w=i->v[i],M,V,cr);
    q=#v+1;
    M=matrix(2*q,2*q);
    for(j=1,2*q,M[1,j]=F(s(j),0));
    for(i=1,q-1,for(j=1,2*q,M[2*i  ,j]=F(s(j),w(i))));
    for(i=1,q-1,for(j=1,2*q,M[2*i+1,j]=derivnum(x=w(i),F(s(j),x))));
    for(j=1,2*q,M[2*q,j]=j+1);
    V=vectorv(2*q);
    V[1]=G3(0);
    for(i=1,q-1,V[2*i  ]=G3(w(i)));
    for(i=1,q-1,V[2*i+1]=derivnum(x=w(i),G3(x)));
    V[2*q]=polcoeff(G3(1/'x)+O('x^3),2);
    a=matsolve(M,V);

    ar3=ceil(a*10^7)/10^7;
    cr = polsturm(numerator(G3('x)-vector(#ar3,j,F(1+j/2,'x))*ar3),0);
    if (cr > 0, warning("F and g are crossing"); return([]););

    ploth(g=0,v[#v],sum(j=1,#ar3,ar3[j]*F(s(j),g))-G3(g));
    sum(j=1,#ar3,ar3[j]*[1.,2.*(1/s(j)+1/(s(j)-1)),psi(s(j)/2),psi((s(j)+1)/2)])
}

A3=H3(5,[0.84,2.04,4.01,9.61])
\\%303 =
\\ [
\\  0.5154369999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999994177,
\\  9.341863310548340548340548340548340548340548340548340548340548340548340548340548340548340548340548340548340548341370,
\\ -1.009487099835328378865408731425389024479455570527936726740252903073170153156903277163802337895364718445415891882005,
\\ -0.2970195609459573108152024241518272992600801464865091978413158315035955095036487378205100871511491143041660944141180
\\ ]
[A3[1], -A3[1]*log(Pi)+max(A3[3],(A3[3]+A3[4])/2),A3[2]]
\\%25 =
\\ [
\\  0.5154369999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999994177,
\\ -1.243289468563200122400271341487974685025114554991852412777051702378767107026401107577587640084244504843591824733118,
\\  9.341863310548340548340548340548340548340548340548340548340548340548340548340548340548340548340548340548340548341370
\\ ]

{
ar3 =
[
    149178011,
  -1773766184,
  11465438478,
 -45115091060,
 114102793523,
-189514259129,
 205612934195,
-140312989024,
  54661946795,
  -9271031235
]~/10^7
};

S(ar,n)=sum(j=1,#ar,ar[j]/n^s(j));
vector(#ar3\2,i,my(k=2*i-1);ceil(solve(x=.0001,10^5,ar3[k]/x^s(k)+ar3[k+1]/x^s(k+1))))
\\%306 = [142, 16, 3, 1, 1]
for(n=2,150,if(S(ar3,n)<0,print(n)))



\\==================================================
\\Lower bound for |B_K| = (1/2)\sum_{\gamma} 4/(1+4\gamma^2) \geq 0.4512\log\disc_\K - 5.2554 n_\K + 5.2784
\\==================================================
\p 115
F(s,g)= 4*(2*s-1)/((2*s-1)^2+4*g^2);
G3(g)=2/(1+4*g^2);
s(j)=1+j/2;
H4(q,v)=
{
    my(w=i->v[i],M,V,cr);
    q=#v+1;
    M=matrix(2*q,2*q);
    for(j=1,2*q,M[1,j]=F(s(j),0)*1.35);
    for(i=1,q-1,for(j=1,2*q,M[2*i  ,j]=F(s(j),w(i))));
    for(i=1,q-1,for(j=1,2*q,M[2*i+1,j]=derivnum(x=w(i),F(s(j),x))));
    for(j=1,2*q,M[2*q,j]=j+1);
    V=vectorv(2*q);
    V[1]=G3(0);
    for(i=1,q-1,V[2*i  ]=G3(w(i)));
    for(i=1,q-1,V[2*i+1]=derivnum(x=w(i),G3(x)));
    V[2*q]=polcoeff(G3(1/'x)+O('x^3),2);
    a=matsolve(M,V);

    ar4=floor(a*10^7)/10^7;
    cr = polsturm(numerator(G3('x)-vector(#ar4,j,F(1+j/2,'x))*ar4),0);
    if (cr > 0, warning("F and g are crossing"); return([]););

    ploth(g=0,v[#v],g^2*(sum(j=1,#ar4,ar4[j]*F(s(j),g))-G3(g)));
    sum(j=1,#ar4,ar4[j]*[1.,2.*(1/s(j)+1/(s(j)-1)),psi(s(j)/2),psi((s(j)+1)/2)])
}

A4=H4(5,[0.84,2.04,4.01,9.61])
\\%133 =
\\ [
\\   0.4512055000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000009,
\\   5.278450823044733044733044733044733044733044733044733044733044733044733044733044733044733044733044733044733044733053,
\\  -0.8057255843033801307073476427948153589596476317855476243632043688304384202671594783728115302191234611163802835346619,
\\  -0.2377168219382190409952743186668983906274048910791531971516430365506601772862371863570633730977376636733211889552825
\\ ]
[A4[1], -A4[1]*log(Pi)+min(A4[3],(A4[3]+A4[4])/2),A4[2]]
\\%134 =
\\ [
\\   0.4512055000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000009,
\\  -1.322234004813001660981819852575747891477821111494407239643794423605560078796927372114621159695410932424919471866219,
\\   5.278450823044733044733044733044733044733044733044733044733044733044733044733044733044733044733044733044733044733053
\\ ]
{
ar4 =
[
   35520843,
 -172090428,
  546191277,
-1118014046,
 1410926499,
 -924631608,
    4105241,
  464624197,
 -311256372,
   69136452
]~/10^7
};
S(ar,n)=sum(j=1,#ar,ar[j]/n^s(j));
my(i=1,k=2*i-1);ceil(solve(x=1,10^5,ar4[k]/x^s(k)+ar4[k+1]/x^s(k+1)))
my(i=2,k=2*i-1);ceil(solve(x=1,10^5,ar4[k]/x^s(k)+ar4[k+1]/x^s(k+1)))
my(i=3,k=2*i-1);ceil(solve(x=0.1,10^5,ar4[k]/x^s(k)+ar4[k+1]/x^s(k+1)))
my(i=4,k=2*i-1);ceil(solve(x=1,10^5,ar4[k]/x^s(k)+ar4[k+1]/x^s(k+1)+ar4[k+2]/x^s(k+2)+ar4[k+3]/x^s(k+3)))     \\  <--- Always positive
\\%306 = [24, 5, 1, "error"]
for(n=1,100,if(S(ar4,n)<0,print(n)))
\\ All S(n) are positive
sum(j=1,#ar4,2*ar4[j]*derivnum(w=s(j),log(zeta(w)))) + (-A4[1]*log(Pi)+min(A4[3],(A4[3]+A4[4])/2))
\\ %592 = -5.255369077550807977060601105886986957553605812948856282344791217341649851345704710317830856863801040




\\ ==================================================
\\Lower bound for |B_K| = (1/2)\sum_{\gamma} 4/(1+4\gamma^2) \geq 0.4840\log\disc_\K - 7.9606 n_\K + 7.9836
\\ ==================================================
\p 100
F(s,g)= (2*s-1)/((s-1/2)^2+g^2);
G(g)=2/(1+4*g^2);
s(j)=1+j/2;
H(v)=
{
    my(q=#v+1,w=i->v[i],s=j->1+j/2,M,V,cr);
    M=matrix(2*q,2*q);
    for(j=1,2*q,M[1,j]=F(s(j),0)*1.1);
    for(i=1,q-1,for(j=1,2*q,M[2*i  ,j]=F(s(j),w(i))));
    for(i=1,q-1,for(j=1,2*q,M[2*i+1,j]=derivnum(x=w(i),F(s(j),x))));
    for(j=1,2*q,M[2*q,j]=j+1);
    V=vector(2*q);
    V[1]=G(0);
    for(i=1,q-1,V[2*i  ]=G(w(i)));
    for(i=1,q-1,V[2*i+1]=derivnum(x=w(i),G(x)));
    V[2*q]=polcoeff(G(1/'x)+O('x^3),2);
    a=matsolve(M,V~);

    ar=floor(a*10^7)/10^7;
    cr = polsturm(numerator(G('x)-vector(#ar,j,F(1+j/2,'x))*ar),0);
    if (cr > 0, warning("F and g are crossing"); return([]););

    ploth(g=0,v[#v],g^2*(sum(j=1,#ar,ar[j]*F(s(j),g))-G(g)));
    sum(j=1,#ar,ar[j]*[1.,2.*(1/s(j)+1/(s(j)-1)),psi(s(j)/2),psi((s(j)+1)/2)])
}

A=H([0.471,1.61,3.,4.99,9.6,13])
\\ 0
\\ %417 = [0.4840060000000000000000000000, 7.983664527713397713397713397713397713397713397713397713397713397713397713397713397713397713397713398, -0.9307393043459929436551252701522079728151965860870365962162906293023996138595964666136367413559279666, -0.2732641517439925202838908792397943056048186983258786828898944456614917523539571465188002936655945735]
[A[1], -A[1]*log(Pi)+min(A[3],(A[3]+A[4])/2),A[2]]
\\ %418 = [0.4840060000000000000000000000, -1.484795437476417724341588968771196507604757159306924888698313277633343127310471598185389105908886638, 7.983664527713397713397713397713397713397713397713397713397713397713397713397713397713397713397713398]

{
    ar ==
[    109617513,
   -1231529756,
    8487137304,
  -39209340826,
  126704121804,
 -291654809645,
  477762626373,
 -542915757926,
  395568630794,
 -135366272755,
  -42005735236,
   69033311406,
  -30214107922,
    4936948932]~/10^7
}

my(i=1,k=2*i-1);ceil(solve(x=1,10^5,ar[k]/x^(1+k/2)+ar[k+1]/x^(1+(k+1)/2)))
my(i=2,k=2*i-1);ceil(solve(x=1,10^5,ar[k]/x^(1+k/2)+ar[k+1]/x^(1+(k+1)/2)))
my(i=3,k=2*i-1);ceil(solve(x=0.1,10^5,ar[k]/x^(1+k/2)+ar[k+1]/x^(1+(k+1)/2)))
my(i=4,k=2*i-1);ceil(solve(x=0.1,10^5,ar[k]/x^(1+k/2)+ar[k+1]/x^(1+(k+1)/2)))
my(i=5,k=2*i-1);ceil(solve(x=0.2,0.3,ar[k]/x^(1+k/2)+ar[k+1]/x^(1+(k+1)/2)+ar[k+2]/x^(1+(k+2)/2)))
my(i=6,k=2*i);ceil(solve(x=0.1,10^5,ar[k]/x^(1+k/2)+ar[k+1]/x^(1+(k+1)/2)))
\\ [ 127, 22, 6, 2, 1, 1]
S(n)=sum(j=1,#ar,ar[j]/n^s(j));
for(n=1,128,if(S(n)<0,print(n)))
\\ All S(n) are positive
sum(j=1,#ar,2*ar[j]*derivnum(w=s(j),log(zeta(w)))) + (-A[1]*log(Pi)+min(A[3],(A[3]+A[4])/2))
\\ %592 = -7.960584432155736532763433679624842597198595180048836227758945735369004600943914932894123948908119260
blrd(x,n=4) = floor(x*10^n)*10.^-n
precision(apply(blrd,[A[1], sum(j=1,#ar,2*ar[j]*derivnum(w=s(j),log(zeta(w)))) + (-A[1]*log(Pi)+min(A[3],(A[3]+A[4])/2)),A[2]]),9)
\\ %421 = [0.4840, -7.9606, 7.9836]