log2_factorial.sage raw

   1  import bisect
   2  
   3  INPUT_BITS = 32
   4  TABLE_BITS = 5
   5  INT_BITS = 64
   6  EXACT_FPBITS = 256
   7  
   8  F = RealField(100) # overkill
   9  
  10  def BestOverApproxInvLog2(mulof, maxd):
  11      """
  12      Compute denominator of an approximation of 1/log(2).
  13  
  14      Specifically, find the value of d (<= maxd, and a multiple of mulof)
  15      such that ceil(d/log(2))/d is the best approximation of 1/log(2).
  16      """
  17      dist=1
  18      best=0
  19      # Precomputed denominators that lead to good approximations of 1/log(2)
  20      for d in [1, 2, 9, 70, 131, 192, 445, 1588, 4319, 11369, 18419, 25469, 287209, 836158, 3057423, 8336111, 21950910, 35565709, 49180508, 161156323, 273132138, 385107953, 882191721]:
  21          kd = lcm(mulof, d)
  22          if kd <= maxd:
  23              n = ceil(kd / log(2))
  24              dis = F((n / kd) - 1 / log(2))
  25              if dis < dist:
  26                  dist = dis
  27                  best = kd
  28      return best
  29  
  30  
  31  LOG2_TABLE = []
  32  A = 0
  33  B = 0
  34  C = 0
  35  D = 0
  36  K = 0
  37  
  38  def Setup(k):
  39      global LOG2_TABLE, A, B, C, D, K
  40      K = k
  41      LOG2_TABLE = []
  42      for i in range(2 ** TABLE_BITS):
  43          LOG2_TABLE.append(int(floor(F(K * log(1 + i / 2**TABLE_BITS, 2)))))
  44  
  45      # Maximum for (2*x+1)*LogK2(x)
  46      max_T = (2^(INPUT_BITS + 1) - 1) * (INPUT_BITS*K - 1)
  47      # Maximum for A
  48      max_A = (2^INT_BITS - 1) // max_T
  49      D = BestOverApproxInvLog2(2 * K, max_A * 2 * K)
  50      A = D // (2 * K)
  51      B = int(ceil(F(D/log(2))))
  52      C = int(floor(F(D*log(2*pi,2)/2)))
  53  
  54  def LogK2(n):
  55      assert(n >= 1 and n < (1 << INPUT_BITS))
  56      bits = Integer(n).nbits()
  57      return K * (bits - 1) + LOG2_TABLE[((n << (INPUT_BITS - bits)) >> (INPUT_BITS - TABLE_BITS - 1)) - 2**TABLE_BITS]
  58  
  59  def Log2Fact(n):
  60      # Use formula (A*(2*x+1)*LogK2(x) - B*x + C) / D
  61      return (A*(2*n+1)*LogK2(n) - B*n + C) // D + (n < 3)
  62  
  63  RES = [int(F(log(factorial(i),2))) for i in range(EXACT_FPBITS * 10)]
  64  
  65  best_worst_ratio = 0
  66  
  67  for K in range(1, 10000):
  68      Setup(K)
  69      assert(LogK2(1) == 0)
  70      assert(LogK2(2) == K)
  71      assert(LogK2(4) == 2 * K)
  72      good = True
  73      worst_ratio = 1
  74      for i in range(1, EXACT_FPBITS * 10):
  75          exact = RES[i]
  76          approx = Log2Fact(i)
  77          if not (approx <= exact and ((approx == exact) or (approx >= EXACT_FPBITS and exact >= EXACT_FPBITS))):
  78              good = False
  79              break
  80          if worst_ratio * exact > approx:
  81              worst_ratio = approx / exact
  82      if good and worst_ratio > best_worst_ratio:
  83          best_worst_ratio = worst_ratio
  84          print("Formula: (%i*(2*x+1)*floor(%i*log2(x)) - %i*x + %i) / %i; log(max_ratio)=%f" % (A, K, B, C, D, RR(-log(worst_ratio))))
  85          print("LOG2K_TABLE: %r" % LOG2_TABLE)
  86