<jupyter_start><jupyter_text>07 Numbers* *Computational Physics*: Ch 2.4, 2.5, 3* Python Tutorial [Floating Point Arithmetic: Issues and Limitations](https://docs.python.org/3/tutorial/floatingpoint.html)  Binary representationComputers store information with two-state logic. This can be represented in a binary system with numbers 0 and 1 (i.e. base 2)Any number can be represented in any base as a polynomial (possibly with infinitely many terms): the digits are $0 \leq x_k < b$ and determine the contribution of the base $b$ raise to the $k$th power.$$q_b = \sum_{k=-\infty}^{+\infty} x_k b^k$$  Integers  Convert 10 (base 10, i.e. $1 \times 10^1 + 0\times 10^0$) into binary (Note: `divmod(x, 2)` is `x // 2, x % 2`, i.e. integer division and remainder):<jupyter_code>divmod(10, 2)
divmod(5, 2)
divmod(2, 2)<jupyter_output><empty_output><jupyter_text>The binary representation of $10_{10}$ is $1010_2$ (keep dividing until there's only 1 left, then collect the 1 and all remainders in reverse order, essentially long division). Double check by multiplying out $1010_2$:<jupyter_code>1*2**3 + 0*2**2 + 1*2**1 + 0*2**0<jupyter_output><empty_output><jupyter_text>or in Python<jupyter_code>int('0b1010', 2)<jupyter_output><empty_output><jupyter_text>Summary: Integers in binary representation**All integers are exactly representable in base 2 with a finite number of digits**.* The sign (+ or –) is represented by a single bit (0 = +, 1 = –). * The number of available "bits" (digits) determines the largest representable integer.  For example, with 8 bits available (a "*byte*"), what is the largest and smallest integer?<jupyter_code>0b1111111  # 7 bits for number, 1 for sign (not included)
-0b1111111<jupyter_output><empty_output><jupyter_text>Sidenote: using numpy to quickly convert integersIf you want to properly sum all terms, use numpy arrays and the element-wise operations:<jupyter_code>import numpy as np
nbits = 7
exponents = np.arange(nbits)
bases = 2*np.ones(nbits)  # base 2
digits = np.ones(nbits)   # all 1, for 1111111 (127 in binary)
exponents, bases, digits
np.sum(digits * bases**exponents)<jupyter_output><empty_output><jupyter_text>Examples: limits of integersWhat is the smallest and largest integer that you can represent1. if you have 4 bits available and only consider non-negative ("unsigned") integers?2. if you have 32 bits and consider positive and negative integers?3. if you have 64 bits and consider positive and negative integers? Smallest and largest 4 bit unsigned integer:<jupyter_code>0b0000
0b1111<jupyter_output><empty_output><jupyter_text>Smallest and largest 32-bit signed integer (int32): 1 bit is sign, 31 bits are available, so the highest number has 31 ones (111...11111). The *next highest* number is 1000...000, a one with 32 bits and 31 zeroes, i.e., $2^{31}$. Thus, the highest number is $2^{31} - 1$:<jupyter_code>2**31 - 1<jupyter_output><empty_output><jupyter_text>(and the smallest number is just $-(2^{31} - 1)$) And int64 (signed):<jupyter_code>max64 = 2**(64-1) - 1
print(-max64, max64)<jupyter_output>-9223372036854775807 9223372036854775807<jupyter_text>Python's arbitrary precision integers  In Python, integers *have arbitrary precision*: integer arithmetic (`+`, `-`, `*`, `//`) is exact and will not overflow. Thus the following code will run forever (until memory is exhausted); if you run it, you can stop the evaluation with the ''Kernel / Interrupt'' menu command in the notebook and then investigate `n` and `nbits`:<jupyter_code>n = 1 
nbits = 1
while True:
    n *= 2
    nbits += 1
type(n)
int.bit_length(n)
nbits<jupyter_output><empty_output><jupyter_text>NumPy has fixed precision integersNumPy data types (dtypes) are fixed precision. Overflows "wrap around":<jupyter_code>import numpy as np
np.array([2**15-1], dtype=np.int16)
np.array([2**15], dtype=np.int16)
np.array([2**15 + 1], dtype=np.int16)<jupyter_output><empty_output><jupyter_text>Binary fractionsDecimal fractions can be represented as binary fractions:Convert $0.125_{10}$ to base 2:<jupyter_code>0.125 * 2  # 0.0
_ * 2      # 0.00
_ * 2      # 0.001<jupyter_output><empty_output><jupyter_text>Thus the binary representation of $0.125_{10}$ is $0.001_2$. General recipe:- multiply by 2- if you get a number < 1, add a digit 0 to the right- if you get a number ≥ 1, add a digit 1 to the right and then use the remainder in the same fashion<jupyter_code>0.3125 * 2    # 0.0
_ * 2         # 0.01
(_ - 1) * 2   # 0.010
_ * 2         # 0.0101<jupyter_output><empty_output><jupyter_text>Thus, 0.3125 is $0.0101_2$. What is the binary representation of decimal $0.1 = \frac{1}{10}$?<jupyter_code>0.1 * 2  # 0.0
_ * 2   # 0.00
_ *  2  # 0.000
_ * 2  # 0.0001
(_ - 1) * 2  # 0.00011
(_ - 1) * 2  # 0.000110
_ * 2  # 0.0001100 
_ * 2  # 0.00011001<jupyter_output><empty_output><jupyter_text>... etc: this is an infinitely repeating fraction and the binary representation of $0.1_{10}$ is $0.000 1100 1100 1100 ..._2$.**Thus, with a finite number of bits, 0.1 is not exactly representable in the computer.** The number 0.1 is not stored exactly in the computer. `print` only shows you a convenient approximation:<jupyter_code>print(0.1)
print("{0:.55f}".format(0.1))<jupyter_output>0.1000000000000000055511151231257827021181583404541015625<jupyter_text>Problems with floating point arithmeticOnly a subset of all real numbers can be represented with **floating point numbers of finite bit size**. Almost all floating point numbers are not exact:<jupyter_code>0.1 + 0.1 + 0.1 == 0.3<jupyter_output><empty_output><jupyter_text>... which should have yielded `True`! But because the machine representation of 0.1 is not exact, the equality cannot be fulfilled.  Representation of floats: IEEE 754Floating point numbers are stored in "scientific notation": e.g. $c = 2.88792458 \times 10^8$ m/s  * **mantissa**: $2.88792458$  * **exponent**: $+8$  * **sign**: +Format: $$x = (-1)^s \times 1.f \times 2^{e - \mathrm{bias}}$$($f$ is $M$ bits long. The leading 1 in the mantissa is assumed and not stored: "ghost" or "phantom" bit.) Format: $$x = (-1)^s \times 1.f \times 2^{e - \mathrm{bias}}$$Note: * In IEEE 754, the highest value of $e$ in the exponent is reserved and not used, e.g. for a 32-bit *float* (see below) the exponent has $(30 - 23) + 1 = 8$ bit and hence the highest number for $e$ is $(2^8 - 1) - 1 = 255 - 1 = 254$. Taking the *bias* into account (for *float*, *bias* = 127), the largest value for the exponent is $2^{254 - 127} = 2^{127}$.* The case of $e=0$ is also special. In this case, the format is $$x = (-1)^s \times 0.f \times 2^{-\mathrm{bias}}$$ i.e. the "ghost 1" becomes a zero, gaining a additional order of magnitude.  IEEE float (32 bit)IEEE *float* uses **32 bits**  * $\mathrm{bias} = 127_{10}$  * bits        sef    bit position3130–2322–0      * **six or seven decimal places of significance** (1 in $2^{23}$)  * range: $1.4 \times 10^{-45} \leq |x_{(32)}| \leq 3.4 \times 10^{38}$<jupyter_code>1/2**23<jupyter_output><empty_output><jupyter_text>IEEE double (64 bit)Python floating point numbers are 64-bit doubles. NumPy has dtypes `float32` and `float64`.IEEE *double* uses **64 bits**  * $\mathrm{bias} = 1023_{10}$  * bits        sef    bit position6362–5251–0      * **about 16 decimal places of significance** (1 in $2^{52}$)  * range: $4.9 \times 10^{-324} \leq |x_{(64)}| \leq 1.8 \times 10^{308}$<jupyter_code>1/2**52<jupyter_output><empty_output><jupyter_text>For numerical calculations, *doubles* are typically required.  Special numbersIEEE 754 also introduces special "numbers" that can result from floating point arithmetic* `NaN` (not a number)* `+INF` and `-INF` (infinity)* `-0` (signed zero) Python itself does not use the IEEE special numbers<jupyter_code>1/0<jupyter_output><empty_output><jupyter_text>But numpy does:<jupyter_code>np.array([1, -1])/np.zeros(2)<jupyter_output>/Users/oliver/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in true_divide
  """Entry point for launching an IPython kernel.<jupyter_text>But beware, you cannot use `INF` to "take limits". It is purely a sign that something bad happened somewhere... And **not a number**, `nan`<jupyter_code>np.zeros(2)/np.zeros(2)<jupyter_output>/Users/oliver/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in true_divide
  """Entry point for launching an IPython kernel.<jupyter_text>Overflow and underflow* underflow: typically just set to zero (and that works well most of the time)* overflow: raises exception or just set to `inf`<jupyter_code>big = 1.79e308
big
2 * big
2 * np.array([big], dtype=np.float64)<jupyter_output>/Users/oliver/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: overflow encountered in multiply
  """Entry point for launching an IPython kernel.<jupyter_text>... but you can (maybe*) just use an even bigger data type:<jupyter_code>2 * np.array([big], dtype=np.float128)<jupyter_output><empty_output><jupyter_text>(*Note: float128 is 80 bit on Linux/macOS and 64 bit on Windows... see [Numpy data types: Extended Precisions](https://docs.scipy.org/doc/numpy/user/basics.types.htmlextended-precision))  Insignificant digits<jupyter_code>x = 1000.2
A = 1000.2 - 1000.0
print(A)
A == 0.2<jupyter_output><empty_output><jupyter_text>... oops<jupyter_code>x = 700
y = 1e-14
x - y
x - y < 700<jupyter_output><empty_output><jupyter_text>... ooops  Machine precisionOnly a limited number of floating point numbers can be represented. This *limited precision* affects calculations:<jupyter_code>x = 5  + 1e-16
x
x == 5<jupyter_output><empty_output><jupyter_text>... oops. **Machine precision** $\epsilon_m$ is defined as the maximum number that can be added to 1 in the computer without changing that number 1:$$1_c + \epsilon_m := 1_c$$Thus, the *floating point representation* $x_c$ of an arbitrary number $x$ is "in the vicinity of $x$"$$x_c = x(1\pm\epsilon), \quad |\epsilon| \leq \epsilon_m$$where we don't know the true value of $\epsilon$. Thus except for powers of 2 (which are represented exactly) **all floating point numbers contain an unknown error in the 6th decimal place (32 bit floats) or 15th decimal (64 bit doubles)**. This error should be treated as a random error because we don't know its magnitude.<jupyter_code>N = 100
eps = 1
for nbits in range(N):
    eps /= 2
    one_plus_eps = 1.0 + eps
    # print("eps = {0}, 1 + eps = {1}".format(eps, one_plus_eps))
    if one_plus_eps == 1.0:
        print("machine precision reached for {0} bits".format(nbits))
        print("eps = {0}, 1 + eps = {1}".format(eps, one_plus_eps))
        break<jupyter_output>machine precision reached for 52 bits
eps = 1.1102230246251565e-16, 1 + eps = 1.0<jupyter_text>Compare to our estimate for the precision of float64:<jupyter_code>1/2**52<jupyter_output><empty_output><jupyter_text>Appendix A quick hack to convert a floating point binary representation to a floating point number.<jupyter_code>bits = "1010.0001100110011001100110011001100110011001100110011"
import math
def bits2number(bits):
    if '.' in bits:
        integer, fraction = bits.split('.')
    else:
        integer = bits
        fraction = ""
    powers = [int(bit) * 2**n for n, bit in enumerate(reversed(integer))]
    powers.extend([int(bit) * 2**(-n) for n, bit in enumerate(fraction, start=1)])
    return math.fsum(powers)
bits2number(bits)
bits2number('1111')
bits2number('0.0001100110011001100110011001100110011001100110011')
bits2number('0.0001100')
bits2number('0.0101')
bits2number("10.10101")
bits2number('0.0111111111111111111111111111111111111111')
bits2number('0.110011001100')<jupyter_output><empty_output><jupyter_text>Python can convert to binary using the `struct` module:<jupyter_code>x = 6.0e-8

import struct
fpack = struct.pack('f', x)    # pack float into bytes
fint = struct.unpack('i', fpack)[0] # unpack to int
m_bits = bin(fint)[-23:]            # mantissa bits
print(m_bits)<jupyter_output>00000001101100101011001<jupyter_text>With phantom bit:<jupyter_code>mantissa_bits = '1.' + m_bits
print(mantissa_bits)
import math
mn, ex = math.frexp(x)
print(mn, ex)<jupyter_output>0.50331648 -23
{'fraction_non_alphanumeric': 0.12651589642740085, 'fraction_numerical': 0.08218616846935431, 'mean_word_length': 4.884763741562199, 'pattern_counts': {'":': 2, '<': 116, '<?xml version=': 0, '>': 113, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 6, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '30', 'dir': 'jupyter-structured-clean-dedup', 'id': '56563', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4834, 'n_tokens_neox': 4079, 'n_words': 1361, 'path': '07_numbers/07-numbers.ipynb', 'repo_name': 'ASU-CompMethodsPhysics-PHY494/PHY494-resources-2019', 'size': '12204'}