(original) (raw)

# Appendix A def dn(n): return lcm([i for i in [1..n]]) def Dmn(m,n): p = factorial(m) for q in prime_range(n+1): while p%q == 0: p = p//q return p # search 1-beta^2 perfect square def search_square(p,e=0): for x in range(1,2^p): beta = x/2^p*2^e if is_square(1-beta^2): print (beta) # search beta^2-1 perfect square def search_square2(p,e=0): for x in range(1,2^p): beta = 1+x/2^p*2^e if is_square(beta^2-1): print (beta) # search beta^2+1 perfect square def search_square3(p,e=0): for x in range(1,2^p): beta = x/2^p*2^e if is_square(beta^2+1): print (beta) # for 2^e <= alpha < 2^(e+1) # KVexp(-3) # [1/8,1/4] Table 4 page 14 # 226 81 41 3.97072324477161978360835994357507e-7714 # KVexp(-2,Kmax=200) # 297 84 59 6.42424857004856442710224968900302e-10107 # KVexp(-1,Kmax=200) # 403 121 59 2.78620107992289202844845627242010e-13733 # KVexp(0,Kmax=200) # 593 170 67 1.62774958625874988306876271573699e-20180 # KVexp(1,Kmax=300) # 920 208 97 8.38784848856714586400121471031909e-31321 # KVexp(2,Kmax=400) # 1485 347 97 5.55969286813273912305565543435633e-50520 def KVexp(e,p=113,Kmax=100): D = 1 A = 2.^max(0,p-1-e) beta_min = exp(RR(2^e)) # we want log(B/2) maximal thus e_beta minimal e_beta = floor(log(beta_min)/log(2.)) B = 2.^max(0,p-1-e_beta) var('E') Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in [1..Kmax]: for L in [2..Lmax]: lhs = K*L*log(E) rhs = D*K*L*log(2.) rhs += D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*1.0) rhs += D*(1+2*log(2.))*(L-1) rhs += D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*1.0) rhs += (K-1)*log(A/2.0) rhs += L*E*2^(e+1.) # L*E*|alpha| rhs += L*1.0*log(E) rhs += (L-1)*log(B/2.0) eq = lhs - rhs # print (K, L, eq) # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1.0) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # convert into ulp(exp(alpha)) y = exp(R(2^(e+1))) maxbound = maxbound/y.ulp() return ceil(-log(maxbound)/log(2.)/p) # return the minimum value of |cos(x)| for x a p-bit floating-point number, # 2^e <= x < 2^(e+1) # CosMin(-1,113) # 0.540302305868139717400936607442977 # CosMin(1,113) # 1.30077151951856715371955660390650e-34 def CosMin(e,p): R = RealField(p) beta_min = abs(cos(R(2^e))) beta_min = min(beta_min, abs(cos(R(2^(e+1)).nextbelow()))) # cos is minimum for x near (1/2+k)*pi kmin = ceil(2^e/pi-1/2) kmax = floor(2^(e+1)/pi-1/2) for k in [kmin..kmax]: x = R((1/2+k)*pi) assert 2^e < x < 2^(e+1), "2^e < x < 2^(e+1)" for j in [-10..10]: xj = x + j*x.ulp() yj = abs(cos(xj)) beta_min = min(beta_min, yj) return beta_min # KVcos(-3,Kmax=200) # 1359 # KVcos(-2,Kmax=400) # 2058 # KVcos(-1,Kmax=400) # 3281 # KVcos(0,Kmax=700,Lmax=200) # 6481 # KVcos(1,Kmax=800) # 10266 # KVcos(2,Kmax=1500,Lmax=300) # 20395 def KVcos(e,p=113,Kmax=100,Lmax=None,step=1,verbose=false): alpha = RIF(2^e,2^(e+1)) # Lemma A.3 D = 2 A = 2.^(2*max(0,p-1-e))*max(1,alpha.abs().upper()) beta = cos(alpha) beta_min = CosMin(e,p) e_beta = floor(log(beta_min)/log(2.)) print ("beta_min=", beta_min, "e_beta=", e_beta) B = 2.^max(0,p-2-e_beta) var('E') if Lmax==None: Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in range(1,Kmax+1,step): for L in range(2,Lmax+1,step): lhs = K*L*log(E) rhs = D*K*L*log(R(2)) rhs += D*(K-1)*log(exp(R(1))*sqrt(R(3)*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*R(1)) rhs += D*(1+2*log(R(2)))*(L-1) rhs += D*log(R(1)*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*R(1)) rhs += (K-1)*log(A/R(2)) rhs += L*E*alpha.abs().upper() # L*E*|alpha'| rhs += L*R(1)*log(E) rhs += (L-1)*log(B/R(2)) eq = lhs - rhs # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # bound1 is sqrt(2*eps) eps1 = maxbound^2/2 # check closest multiple of pi in range kmin = ceil(alpha.lower()/RR(pi)) kmax = floor(alpha.upper()/RR(pi)) if kmin > kmax: delta = 1-beta.abs().upper() else: delta = 1 for k in [kmin..kmax]: x = R(k*pi) for j in [-1..1]: xj = x+j*x.ulp() Xj = xj.exact_rational() y = numerical_approx(cos(Xj),3*p+10) d = 1 - abs(y) if d < delta: delta = d print ("delta=", delta) # bound2 is 2*eps/sqrt(delta) eps2 = maxbound*sqrt(delta)/2 if eps2 > eps1: # print ("Using bound 2*epsilon/sqrt(delta)") eps = eps2 else: eps = eps1 # convert into ulp(cos(alpha)) RIFp = RealIntervalField(p) y = cos(RIFp(alpha)) maxbound = eps/y.abs().upper().ulp() return ceil(-log(maxbound)/log(2.)/p) # return the minimum value of |sin(x)| for x a p-bit floating-point number, # 2^e <= x < 2^(e+1) def SinMin(e,p): R = RealField(p) beta_min = abs(sin(R(2^e))) beta_min = min(beta_min, abs(sin(R(2^(e+1)).nextbelow()))) # sin is minimum for x near k*pi kmin = ceil(2^e/pi) kmax = floor(2^(e+1)/pi) for k in [kmin..kmax]: x = R(k*pi) for j in [-10..10]: xj = x + j*x.ulp() yj = abs(sin(xj)) beta_min = min(beta_min, yj) return beta_min # KVsin(-3,Kmax=200) # 1371 # KVsin(-2,Kmax=200,Lmax=150) # 2070 # KVsin(-1,Kmax=400) # 3288 # KVsin(0,Kmax=600,Lmax=200) # 5678 # KVsin(1,Kmax=800,Lmax=300) # 11408 # KVsin(2,Kmax=1500,Lmax=300) # 20395 def KVsin(e,p=113,Kmin=1,Kmax=100,Lmax=None): alpha = RIF(2^e,2^(e+1)) # Lemma A.3 D = 2 A = 2.^(2*max(0,p-1-e))*max(1,alpha.abs().upper()) beta = sin(alpha) beta_min = SinMin(e,p) e_beta = floor(log(beta_min)/log(2.)) print ("beta_min=", beta_min, "e_beta=", e_beta) B = 2.^max(0,p-2-e_beta) var('E') if Lmax==None: Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in [Kmin..Kmax]: for L in [2..Lmax]: lhs = K*L*log(E) rhs = D*K*L*log(2.) rhs += D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*1.0) rhs += D*(1+2*log(2.))*(L-1) rhs += D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*1.0) rhs += (K-1)*log(A/2.0) rhs += L*E*alpha.abs().upper() # L*E*|alpha'| rhs += L*1.0*log(E) rhs += (L-1)*log(B/2.0) eq = lhs - rhs # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1.0) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # bound1 is sqrt(2*eps) eps1 = maxbound^2/2 # check closest odd multiple of pi/2 in range kmin = ceil(alpha.lower()/RR(pi/2)) kmax = floor(alpha.upper()/RR(pi/2)) if is_even(kmin): kmin += 1 if is_even(kmax): kmax -= 1 if kmin > kmax: delta = beta.abs().upper() else: delta = 1 for k in range(kmin,kmax+1,2): x = R(k*pi/2) for j in [-1..1]: xj = x+j*x.ulp() Xj = xj.exact_rational() y = numerical_approx(sin(Xj),3*p+10) d = 1 - abs(y) if d < delta: delta = d print ("delta=", delta) # bound2 is 2*eps/sqrt(delta) eps2 = maxbound*sqrt(delta)/2 if eps2 > eps1: print ("Using bound 2*epsilon/sqrt(delta)") eps = eps2 else: eps = eps1 # convert into ulp(sin(alpha)) Rp = RealIntervalField(p) y = sin(Rp(alpha)) maxbound = eps/y.abs().upper().ulp() return ceil(-log(maxbound)/log(2.)/p) # KVcosh(-3,Kmax=200) # 682 # KVcosh(-2,Kmax=300) # 1057 # KVcosh(-1,Kmax=400) # 1695 # KVcosh(0,Kmax=500) # 2889 # KVcosh(1,Kmax=900) # 5285 # KVcosh(2,Kmax=1700,Lmax=150) # 10155 def KVcosh(e,p=113,Kmax=100,Lmax=None,step=1): alpha = RIF(2^e,2^(e+1)) # Lemma A.4 D = 2 A = 2.^(max(0,p-1-e))*max(1,alpha.abs().upper()) beta = cosh(alpha) beta_min = beta.abs().lower() e_beta = floor(log(beta_min)/log(2.)) B = 2.^max(0,p-2-e_beta)*2*beta.abs().upper() var('E') if Lmax==None: Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in range(1,Kmax+1,step): for L in range(2,Lmax+1,step): lhs = K*L*log(E) rhs = D*K*L*log(2.) rhs += D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*1.0) rhs += D*(1+2*log(2.))*(L-1) rhs += D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*1.0) rhs += (K-1)*log(A/2.0) rhs += L*E*alpha.abs().upper() # L*E*|alpha'| rhs += L*1.0*log(E) rhs += (L-1)*log(B/2.0) eq = lhs - rhs # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1.0) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # bound1 is sqrt(2*eps) eps1 = maxbound^2/2 # bound2 is 2*eps/sqrt(delta) delta = beta.abs().lower() - 1 eps2 = maxbound*sqrt(delta)/2 if eps2 > eps1: print ("Using bound 2*epsilon/sqrt(delta)") eps = eps2 else: eps = eps1 # convert into ulp(cosh(alpha)) Rp = RealIntervalField(p) y = cosh(Rp(alpha)) maxbound = eps/y.abs().upper().ulp() return ceil(-log(maxbound)/log(2.)/p) # KVsinh(-3,Kmax=200) # 688 # KVsinh(-2,Kmax=300,Lmax=100) # 1062 # KVsinh(-1,Kmax=400,Lmax=150) # 1698 # KVsinh(0,Kmax=600,Lmax=200) # 2889 # KVsinh(1,Kmax=900,Lmax=200) # 5285 # KVsinh(2,Kmax=1700,Lmax=150) # 10155 def KVsinh(e,p=113,Kmax=100,Lmax=None,step=1): alpha = RIF(2^e,2^(e+1)) # Lemma A.5 D = 2 A = 2.^max(0,p-1-e)*max(1,alpha.abs().upper()) beta = sinh(alpha) beta_min = beta.abs().lower() beta_max = beta.abs().upper() e_beta = floor(log(beta_min)/log(2.)) B = 2.^max(0,p-2-e_beta)*(beta_max + sqrt(beta_max^2+1)) var('E') if Lmax==None: Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in range(1,Kmax+1,step): for L in range(2,Lmax+1,step): lhs = K*L*log(E) rhs = D*K*L*log(2.) rhs += D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*1.0) rhs += D*(1+2*log(2.))*(L-1) rhs += D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*1.0) rhs += (K-1)*log(A/2.0) rhs += L*E*alpha.abs().upper() # L*E*|alpha'| rhs += L*1.0*log(E) rhs += (L-1)*log(B/2.0) eq = lhs - rhs # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1.0) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # bound is 4*eps eps = maxbound/4 # convert into ulp(sinh(alpha)) Rp = RealIntervalField(p) y = sinh(Rp(alpha)) maxbound = eps/y.abs().upper().ulp() return ceil(-log(maxbound)/log(2.)/p) # return the minimum value of |tan(x)| for x a p-bit floating-point number, # 2^e <= x < 2^(e+1) def TanMin(e,p): R = RealField(p) beta_min = abs(tan(R(2^e))) beta_min = min(beta_min, abs(tan(R(2^(e+1)).nextbelow()))) # tan is minimum for x near k*pi kmin = ceil(2^e/pi) kmax = floor(2^(e+1)/pi) for k in [kmin..kmax]: x = R(k*pi) for j in [-10..10]: xj = x + j*x.ulp() yj = abs(tan(xj)) beta_min = min(beta_min, yj) return beta_min # return the maximum value of |tan(x)| for x a p-bit floating-point number, # 2^e <= x < 2^(e+1) # TanMax(0,113,verbose=true) # x= 0x1.921fb54442d18469898cc51701b8p+0 y= 0x1.1c46bd57277993a2ee60193c957bp+114 # 2.30632355873715617276619838163737e34 def TanMax(e,p,verbose=false): R = RealField(p) R2 = RealField(p+10) beta_max = abs(tan(R(2^e))) beta_max = max(beta_max, abs(tan(R(2^(e+1)).nextbelow()))) # tan is maximum for x near (1/2+k)*pi kmin = ceil(2^e/pi-1/2) kmax = floor(2^(e+1)/pi-1/2) for k in [kmin..kmax]: x = R((1/2+k)*pi) assert 2^e < x < 2^(e+1), "2^e < x < 2^(e+1)" for j in [-10..10]: xj = x + j*x.ulp() yj = abs(tan(xj)) if yj > beta_max: beta_max = yj if verbose: # print ("x=", get_hex(xj), "y=", get_hex(yj)) print ("x=", R2(xj), "y=", R2(yj)) return beta_max # KVtan(-3,Kmax=200) # 303 # KVtan(-2,Kmax=300,Lmax=100) # 410 # KVtan(-1,Kmax=300,Lmax=100) # 604 # KVtan(0,Kmax=400,Lmax=150) # 1194 # KVtan(1,Kmax=500,Lmax=200) # 1854 # KVtan(2,Kmax=800,Lmax=200) # 3361 769 97 2.07051030459971451267203231705123e-114319 # GP: 3359 E=34.0081009380320748490947240206690 def KVtan(e,p=113,Kmax=100,Lmax=None,step=1): alpha = RIF(2^e,2^(e+1)) # Lemma A.6 D = 1 A = 2.^max(0,p-1-e) beta = tan(alpha) beta_min = TanMin(e,p) print ("beta_min=", beta_min) beta_max = TanMax(e,p) print ("beta_max=", beta_max) e_beta = floor(log(beta_min)/log(2.)) B = 2.^max(0,p-1-e_beta)*sqrt(beta_max^2+1) print("log(A) = ",log(A)) print("log(B) = ",log(B)) var('E') if Lmax==None: Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in range(1,Kmax+1,step): for L in range(2,Lmax+1,step): lhs = K*L*log(E) t1 = D*K*L*log(2.) rhs = t1 t2 = D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += t2 t3 = D*log(Dmn(K-1,L-1)*1.0) rhs += t3 t4 = D*(1+2*log(2.))*(L-1) rhs += t4 t5 = D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += t5 t6 = log(factorial(K-1)*1.0) rhs += t6 t7 = (K-1)*log(A/2.0) rhs += t7 t8 = (L-1)*log(B/2.0) rhs += t8 # temoin = rhs rhs += L*E*2*alpha.abs().upper() # |alpha'| = 2*|alpha| rhs += L*1.0*log(E) eq = lhs - rhs # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1.0) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: # print(temoin, E0) print (K, L, bound) maxbound = bound if maxbound==0: return infinity # bound is 2*eps eps = maxbound/2 # convert into ulp(tan(alpha)) Rp = RealIntervalField(p) maxbound = eps/R(beta_max).ulp() return ceil(-log(maxbound)/log(2.)/p) # return the maximum value of |cot(x)| for x a p-bit floating-point number, # 2^e <= x < 2^(e+1) # sage: CotMax(1,113) # 1.15316177936857808638309919081869e34 def CotMax(e,p): return 1/TanMin(e,p) # return the minimum value of |cot(x)| for x a p-bit floating-point number, # 2^e <= x < 2^(e+1) # sage: CotMin(1,113) # 0.457657554360285763750277410432047 def CotMin(e,p,verbose=false): return 1/TanMax(e,p) # KVcot(-3,Kmax=200) # 302 # KVcot(-2,Kmax=300,Lmax=100) # 410 # KVcot(-1,Kmax=300,Lmax=100) # 604 # KVcot(0,Kmax=400,Lmax=100) # 1196 # KVcot(1,Kmax=500,Lmax=200) # 1855 417 97 1.23631646613908257026933496750760e-63079 # KVcot(2,Kmax=900,Lmax=200) # 3360 769 97 4.78383716255412934762708848833450e-114278 def KVcot(e,p=113,Kmax=100,Lmax=None,step=1): alpha = RIF(2^e,2^(e+1)) # Lemma A.6 D = 1 A = 2.^max(0,p-1-e) beta = cot(alpha) beta_min = CotMin(e,p) print ("beta_min=", beta_min) beta_max = CotMax(e,p) print ("beta_max=", beta_max) e_beta = floor(log(beta_min)/log(2.)) B = 2.^max(0,p-1-e_beta)*sqrt(beta_max^2+1) var('E') if Lmax==None: Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in range(1,Kmax+1,step): for L in range(2,Lmax+1,step): lhs = K*L*log(E) rhs = D*K*L*log(2.) rhs += D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*1.0) rhs += D*(1+2*log(2.))*(L-1) rhs += D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*1.0) rhs += (K-1)*log(A/2.0) rhs += L*E*2*alpha.abs().upper() # |alpha'| = 2*|alpha| rhs += L*1.0*log(E) rhs += (L-1)*log(B/2.0) eq = lhs - rhs # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1.0) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # bound is 2*eps eps = maxbound/2 # convert into ulp(cot(alpha)) Rp = RealIntervalField(p) maxbound = eps/beta_max.ulp() return ceil(-log(maxbound)/log(2.)/p) # KVtanh(-3,Kmax=200) # 303 # KVtanh(-2,Kmax=200) # 409 # KVtanh(-1,Kmax=300) # 600 # KVtanh(0,Kmax=300) # 931 # KVtanh(1,Kmax=400,Lmax=150) # 1507 # KVtanh(2,Kmax=700) # 2634 def KVtanh(e,p=113,Kmax=100,Lmax=None,step=1): alpha = RIF(2^e,2^(e+1)) # Lemma A.7 D = 1 A = 2.^max(0,p-2-e) beta = tanh(alpha) beta_min = beta.abs().lower() beta_max = beta.abs().upper() e_beta = floor(log(beta_min)/log(2.)) B = 2.^max(0,p-1-e_beta)*(beta_max+1) var('E') if Lmax==None: Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in range(1,Kmax+1,step): for L in range(2,Lmax+1,step): lhs = K*L*log(E) rhs = D*K*L*log(2.) rhs += D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*1.0) rhs += D*(1+2*log(2.))*(L-1) rhs += D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*1.0) rhs += (K-1)*log(A/2.0) rhs += L*E*2*alpha.abs().upper() # |alpha'| = 2*|alpha| rhs += L*1.0*log(E) rhs += (L-1)*log(B/2.0) eq = lhs - rhs # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1.0) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # bound is 2*eps eps = maxbound/2 # convert into ulp(tanh(alpha)) Rp = RealIntervalField(p) y = tanh(Rp(alpha)) maxbound = eps/y.abs().upper().ulp() return ceil(-log(maxbound)/log(2.)/p) # KVcoth(-3,Kmax=200) # 298 # KVcoth(-2,Kmax=200) # 405 # KVcoth(-1,Kmax=200) # 598 # KVcoth(0,Kmax=300) # 929 # KVcoth(1,Kmax=400) # 1504 # KVcoth(2,Kmax=700,Lmax=150) # 2631 def KVcoth(e,p=113,Kmax=100,Lmax=None,step=1): alpha = RIF(2^e,2^(e+1)) # Lemma A.7 D = 1 A = 2.^max(0,p-2-e) # FIXME: max(-1,p-2-e) ? beta = coth(alpha) beta_min = beta.abs().lower() beta_max = beta.abs().upper() e_beta = floor(log(beta_min)/log(2.)) B = 2.^max(0,p-1-e_beta)*(beta_max+1) var('E') if Lmax==None: Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in range(1,Kmax+1,step): for L in range(2,Lmax+1,step): lhs = K*L*log(E) rhs = D*K*L*log(2.) rhs += D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*1.0) rhs += D*(1+2*log(2.))*(L-1) rhs += D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*1.0) rhs += (K-1)*log(A/2.0) rhs += L*E*2*alpha.abs().upper() # |alpha'| = 2*|alpha| rhs += L*1.0*log(E) rhs += (L-1)*log(B/2.0) eq = lhs - rhs # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1.0) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if K in [630,631] and L==97: print (K, L, bound) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # bound is 2*eps eps = maxbound/2 # convert into ulp(coth(alpha)) Rp = RealIntervalField(p) y = coth(Rp(alpha)) maxbound = eps/y.abs().upper().ulp() return ceil(-log(maxbound)/log(2.)/p) # Table 4 # KVlog(-3,Kmax=200) # 609 193 59 4.27418289148934e-20736 # KVlog(-2,Kmax=200) # 480 148 59 9.10218373904143e-16355 # KVlog(-1,Kmax=200) # 340 99 59 3.19010515460611e-11569 # KVlog(0,Kmax=200) # 338 99 59 6.85661671757894e-11516 # KVlog(1,Kmax=200) # 475 146 59 3.11123529765040e-16171 # KVlog(2,Kmax=200) # 599 172 67 1.13109261763636e-20406 def KVlog(e,p=113,Kmax=100): alpha = RIF(2^e,2^(e+1)) alpha_p = log(alpha) D = 1 alpha_p_max = alpha_p.abs().upper() e_alpha_p = floor(log(alpha_p_max)/log(2.)) A = 2.^max(0,p-1-e_alpha_p) beta_p = alpha B = 2.^max(0,p-1-e) var('E') Lmax = Kmax//2 maxbound = 0 for K in [1..Kmax]: for L in [2..Lmax]: lhs = K*L*log(E) rhs = D*K*L*log(2.) rhs += D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*1.0) rhs += D*(1+2*log(2.))*(L-1) rhs += D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*1.0) rhs += (K-1)*log(A/2.0) rhs += L*E*alpha_p_max # L*E*|alpha'| rhs += L*1.0*log(E) rhs += (L-1)*log(B/2.0) eq = lhs - rhs # print (K, L, eq) # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = RR(1.0) E1 = RR(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # apply Lemma 4.8, derivative of exp is itself maxbound = maxbound/exp(alpha_p).abs().upper() # convert into ulp(log(alpha)) R113 = RealIntervalField(p) y = log(R113(2^e,2^(e+1))) maxbound = maxbound/y.abs().upper().ulp() return ceil(-log(maxbound)/log(2.)/p) # KVatan(-3,Kmax=200) # 300 84 59 1.11921270606808760583582490725580e-10180 # KVatan(-2,Kmax=200) # 393 117 59 4.11007219629719507198624396843459e-13343 # KVatan(-1,Kmax=200) # 516 160 59 3.03403256326807970573801568957308e-17549 # KVatan(0,Kmax=200) # 631 181 67 2.41324650416169506570688895428861e-21440 # KVatan(1,Kmax=300) # 707 206 67 1.13292993227885438343273724619625e-24040 # KVatan(2,Kmax=250) # 748 209 71 8.45351622049320602544323908775202e-25434 def KVatan(e,p=113,Kmax=100,step=1): alpha = RIF(2^e,2^(e+1)) alpha_p = atan(alpha) alpha_p_max = alpha_p.abs().upper() e_alpha_p = floor(log(alpha_p_max)/log(2.)) # Lemma A.6 D = 1 A = 2.^max(0,p-2-e_alpha_p) beta = alpha beta_max = beta.abs().upper() B = 2.^max(0,p-1-e)*sqrt(beta_max^2+1) var('E') Lmax = Kmax//2 maxbound = 0 R = RealField(p) for K in range(1,Kmax+1,step): for L in range(2,Lmax+1,step): lhs = K*L*log(E) rhs = D*K*L*log(2.) rhs += D*(K-1)*log(exp(1.)*sqrt(3.*L)*dn(L-1)) rhs += D*log(Dmn(K-1,L-1)*1.0) rhs += D*(1+2*log(2.))*(L-1) rhs += D*log(1.0*min(dn(L-2)^(K-1),factorial(L-2))) rhs += log(factorial(K-1)*1.0) rhs += (K-1)*log(A/2.0) rhs += L*E*2*alpha_p.abs().upper() # |alpha'| = 2*|alpha| rhs += L*1.0*log(E) rhs += (L-1)*log(B/2.0) eq = lhs - rhs # eq = a*E + b*log(E) + c a = eq.coefficient(E,1) b = eq.coefficient(log(E),1) c = eq.coefficient(E,0)(E=1) assert eq == a*E + b*log(E) + c, "eq == a*E + b*log(E) + c" assert a < 0, "a < 0" assert b >= 0, "b >= 0" if b==0: # derivative is a E0 = -c/a if E0 <= 1: continue else: # derivative is a + b/E E1 = -b/a # E1 > 0 # check value in E1 is >= 0 eq1 = eq(E=E1) if E1 <= 1 or eq1 < 0: continue # now look for smallest E such that eq = 0 E0 = R(1.0) E1 = R(E1) while E0.nextabove() != E1: c = (E0+E1)/2 if eq(E=c) < 0: E0 = c else: E1 = c E0 = E1 assert E0 > 1, "E0 > 1" bound = E0^(-K*L) if bound > maxbound: print (K, L, bound) maxbound = bound if maxbound==0: return infinity # bound is 2*eps eps = maxbound/2 # apply Lemma 4.8, derivative of tan(x) is 1 + tan(x)^2 dd = tan(alpha_p)^2+1 maxbound = maxbound/dd.abs().upper() # convert into ulp(atan(alpha)) Rp = RealIntervalField(p) y = atan(Rp(alpha)) maxbound = eps/y.abs().upper().ulp() return ceil(-log(maxbound)/log(2.)/p)