compute-powers.py 1.6 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253
  1. #!/usr/bin/env python
  2. # Compute 10 ** exp with exp in the range [min_exponent, max_exponent] and print
  3. # normalized (with most-significant bit equal to 1) significands in hexadecimal.
  4. from __future__ import print_function
  5. min_exponent = -348
  6. max_exponent = 340
  7. step = 8
  8. significand_size = 64
  9. exp_offset = 2000
  10. class fp:
  11. pass
  12. powers = []
  13. for i, exp in enumerate(range(min_exponent, max_exponent + 1, step)):
  14. result = fp()
  15. n = 10 ** exp if exp >= 0 else 2 ** exp_offset / 10 ** -exp
  16. k = significand_size + 1
  17. # Convert to binary and round.
  18. binary = '{:b}'.format(n)
  19. result.f = (int('{:0<{}}'.format(binary[:k], k), 2) + 1) / 2
  20. result.e = len(binary) - (exp_offset if exp < 0 else 0) - significand_size
  21. powers.append(result)
  22. # Sanity check.
  23. exp_offset10 = 400
  24. actual = result.f * 10 ** exp_offset10
  25. if result.e > 0:
  26. actual *= 2 ** result.e
  27. else:
  28. for j in range(-result.e):
  29. actual /= 2
  30. expected = 10 ** (exp_offset10 + exp)
  31. precision = len('{}'.format(expected)) - len('{}'.format(actual - expected))
  32. if precision < 19:
  33. print('low precision:', precision)
  34. exit(1)
  35. print('Significands:', end='')
  36. for i, fp in enumerate(powers):
  37. if i % 3 == 0:
  38. print(end='\n ')
  39. print(' {:0<#16x}'.format(fp.f, ), end=',')
  40. print('\n\nExponents:', end='')
  41. for i, fp in enumerate(powers):
  42. if i % 11 == 0:
  43. print(end='\n ')
  44. print(' {:5}'.format(fp.e), end=',')
  45. print('\n\nMax exponent difference:',
  46. max([x.e - powers[i - 1].e for i, x in enumerate(powers)][1:]))