C4 =
  (X)->d12psi(p(X))/(2038431744*(Pi^8))+11*d8psi(p(X))/(5898240*(Pi^6))+19*d4psi(p(X))/(24576*(Pi^4))+Psi(p(X))/(128*(Pi^2))

N =
  (X)->floor(u(X)^2)

Psi =
  (Z)->cos(2*Pi*(Z^2-Z-1.0/16))/cos(2*Pi*Z)

R1 =
  (X)->(-1)^(N(X)-1)*Psi(p(X))/u(X)

R2 =
  (X)->(-1)^(N(X)-1)*(Psi(p(X))/u(X)-(1/(96*Pi^2))*d3psi(p(X))/(u(X)^3))

R3 =
  (X)->(-1)^(N(X)-1)*(Psi(p(X))/u(X)-(1/(96*Pi^2))*d3psi(p(X))/(u(X)^3)+(d2psi(p(X))/(64*(Pi^2))+d6psi(p(X))/(18432*(Pi^4)))/(u(X)^5))

R4 =
  (X)->(-1)^(N(X)-1)*(Psi(p(X))/u(X)-(1/(96*Pi^2))*d3psi(p(X))/(u(X)^3)+(d2psi(p(X))/(64*(Pi^2))+d6psi(p(X))/(18432*(Pi^4)))/(u(X)^5)-(d1psi(p(X))/(64*(Pi^2))+d5psi(p(X))/(3840*(Pi^4))+d9psi(p(X))/(5308416*(Pi^6)))/(u(X)^7))

R5 =
  (X)->(-1)^(N(X)-1)*(Psi(p(X))/u(X)-(1/(96*Pi^2))*d3psi(p(X))/(u(X)^3)+(d2psi(p(X))/(64*(Pi^2))+d6psi(p(X))/(18432*(Pi^4)))/(u(X)^5)-(d1psi(p(X))/(64*(Pi^2))+d5psi(p(X))/(3840*(Pi^4))+d9psi(p(X))/(5308416*(Pi^6)))/(u(X)^7)+(d12psi(p(X))/(2038431744*(Pi^8))+11*d8psi(p(X))/(5898240*(Pi^6))+19*d4psi(p(X))/(24576*(Pi^4))+Psi(p(X))/(128*(Pi^2)))/(u(X)^9))

Z =
  (X)->2*sum(Y=1,floor(sqrt(X/(2*Pi))),cos(rstheta(X)-X*log(Y))/sqrt(Y))

d12psi =
  (Z)->(Psi(Z-6*epsilon)-12*Psi(Z-5*epsilon)+66*Psi(Z-4*epsilon)-220*Psi(Z-3*epsilon)+495*Psi(Z-2*epsilon)-792*Psi(Z-epsilon)+924*Psi(Z)-792*Psi(Z+epsilon)+495*Psi(Z+2*epsilon)-220*Psi(Z+3*epsilon)+66*Psi(Z+4*epsilon)-12*Psi(Z+5*epsilon)+Psi(Z+6*epsilon))/(epsilon^12)

d1psi =
  (Z)->(-0.5*Psi(Z-epsilon)+0.5*Psi(Z+epsilon))/epsilon

d2psi =
  (Z)->(Psi(Z-epsilon)-2*Psi(Z)+Psi(Z+epsilon))/(epsilon^2)

d3psi =
  (Z)->(-0.5*Psi(Z-epsilon*2)+Psi(Z-epsilon)-Psi(Z+epsilon)+0.5*Psi(Z+epsilon*2))/(epsilon^3)

d4psi =
  (Z)->(Psi(Z-2*epsilon)-4*Psi(Z-epsilon)+6*Psi(Z)-4*Psi(Z+epsilon)+Psi(Z+2*epsilon))/(epsilon^4)

d5psi =
  (Z)->(-0.5*Psi(Z-3*epsilon)+2*Psi(Z-2*epsilon)-2.5*Psi(Z-epsilon)+2.5*Psi(Z+epsilon)-2*Psi(Z+2*epsilon)+0.5*Psi(Z+3*epsilon))/(epsilon^5)

d6psi =
  (Z)->(Psi(Z-3*epsilon)-6*Psi(Z-2*epsilon)+15*Psi(Z-epsilon)-20*Psi(Z)+15*Psi(Z+epsilon)-6*Psi(Z+2*epsilon)+Psi(Z+3*epsilon))/(epsilon^6)

d8psi =
  (Z)->(Psi(Z-4*epsilon)-8*Psi(Z-3*epsilon)+28*Psi(Z-2*epsilon)-56*Psi(Z-epsilon)+70*Psi(Z)-56*Psi(Z+epsilon)+28*Psi(Z+2*epsilon)-8*Psi(Z+3*epsilon)+Psi(Z+4*epsilon))/(epsilon^8)

d9psi =
  (Z)->(-Psi(Z-5*epsilon)+8*Psi(Z-4*epsilon)-27*Psi(Z-3*epsilon)+48*Psi(Z-2*epsilon)-42*Psi(Z-epsilon)+42*Psi(Z+epsilon)-48*Psi(Z+2*epsilon)+27*Psi(Z+3*epsilon)-8*Psi(Z+4*epsilon)+Psi(Z+5*epsilon))/(2*(epsilon^9))

p =
  (X)->u(X)^2-N(X)

rs1 =
  (X)->Z(X)+R1(X)

rs2 =
  (X)->Z(X)+R2(X)

rs3 =
  (X)->Z(X)+R3(X)

rs4 =
  (X)->Z(X)+R4(X)

rs5 =
  (X)->Z(X)+R5(X)

rstheta =
  (X)->(X/2.0)*log(X/(2.0*Pi))-X/2.0-Pi/8.0+1/(48.0*X)+7.0/(5760.0*X^3)+31.0/(80640*X^5)+381.0/(1290240*X^7)

u =
  (X)->(X/(2*Pi))^(0.25)