This repository has been archived by the owner on Aug 19, 2024. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 2
/
ulp.py
165 lines (116 loc) · 5.27 KB
/
ulp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
"""
Documentation: http://neerc.ifmo.ru/wiki/index.php?title=Представление_вещественных_чисел&oldid=84685
"""
import numpy as np
class FloatingTypeError(TypeError):
message: str = "The given floating point type must be either np.floating."
def __init__(self, message: str = message):
super().__init__(message)
def get_machine_epsilon(number_type: np.floating) -> np.floating:
"""Get the machine epsilon for a given floating point type.
Args:
number_type (np.floating): The floating point type to get the machine epsilon for.
Returns:
np.floating: The machine epsilon for the given floating point type.
Raises:
FloatingTypeError: The given floating point type must be either np.floating.
Doctests:
>>> get_machine_epsilon(np.float32) == np.finfo(np.float32).eps
True
>>> get_machine_epsilon(np.float64) == np.finfo(np.float64).eps
True
"""
if not isinstance(number_type(1.0), np.floating):
raise FloatingTypeError()
machine_epsilon = number_type(1.0)
while number_type(1.0) + machine_epsilon > number_type(1.0):
machine_epsilon /= number_type(2.0)
return machine_epsilon * number_type(2.0)
def get_number_of_mantissa_bits(number_type: np.floating) -> int:
"""Get the number of mantissa bits for a given floating point type.
Args:
number_type (np.floating): The floating point type to get the number of mantissa bits for.
Returns:
int: The number of mantissa bits for the given floating point type.
Raises:
FloatingTypeError: The given floating point type must be either np.floating.
Doctests:
>>> get_number_of_mantissa_bits(np.float32) == np.finfo(np.float32).nmant
True
>>> get_number_of_mantissa_bits(np.float64) == np.finfo(np.float64).nmant
True
"""
if not isinstance(number_type(1.0), np.floating):
raise FloatingTypeError()
number_of_mantissa_bits = 0
while number_type(1.0) + number_type(2.0 ** -(number_of_mantissa_bits + 1)) > number_type(1.0):
number_of_mantissa_bits += 1
return number_of_mantissa_bits
def get_maximum_exponent(number_type: np.floating) -> int:
"""Get the maximum exponent for a given floating point type.
Args:
number_type (np.floating): The floating point type to get the maximum exponent for.
Returns:
int: The maximum exponent for the given floating point type.
Raises:
FloatingTypeError: The given floating point type must be either np.floating.
Doctests:
>>> get_maximum_exponent(np.float32) == np.finfo(np.float32).maxexp
True
>>> get_maximum_exponent(np.float64) == np.finfo(np.float64).maxexp
True
"""
if not isinstance(number_type(1.0), np.floating):
raise FloatingTypeError()
maximum_exponent = 0
maximum_exponent_value = number_type(1.0)
while maximum_exponent_value != number_type("inf"):
maximum_exponent_value *= number_type(2.0)
maximum_exponent += 1
return maximum_exponent
def get_minimum_exponent(number_type: np.floating) -> int:
"""Get the minimum exponent for a given floating point type.
Args:
number_type (np.floating): The floating point type to get the minimum exponent for.
Returns:
int: The minimum exponent for the given floating point type.
Raises:
FloatingTypeError: The given floating point type must be either np.floating.
Doctests:
>>> get_minimum_exponent(np.float32) == np.finfo(np.float32).minexp
True
>>> get_minimum_exponent(np.float64) == np.finfo(np.float64).minexp
True
"""
if not isinstance(number_type(1.0), np.floating):
raise FloatingTypeError()
minimum_exponent = 0
minimum_exponent_value = number_type(1.0)
while minimum_exponent_value != number_type(0.0):
minimum_exponent_value /= number_type(2.0)
minimum_exponent -= 1
return minimum_exponent + get_number_of_mantissa_bits(number_type) + 1
if __name__ == "__main__":
for np_type in [np.float32, np.float64]:
epsilon = get_machine_epsilon(np_type)
mantissa_bits = get_number_of_mantissa_bits(np_type)
max_exp = get_maximum_exponent(np_type)
min_exp = get_minimum_exponent(np_type)
print(f"{np_type.__name__}:")
print(f"Machine epsilon: {epsilon}")
print(f"Number of mantissa bits: {mantissa_bits}")
print(f"Maximum exponent: {max_exp}")
print(f"Minimum exponent: {min_exp}")
print()
np_type = np.float32
epsilon = get_machine_epsilon(np_type)
print(f"1 < 1 + eps: {np_type(1.0) < np_type(1.0) + epsilon}")
print(f"1 == 1 + eps / 2: {np_type(1.0) == np_type(1.0) + epsilon / np_type(2.0)}")
print(f"1 < 1 + eps + eps / 2: {np_type(1.0) < np_type(1.0) + epsilon + epsilon / np_type(2.0)}")
print(f"1 + eps/2 < 1 + eps: {np_type(1.0) + epsilon / np_type(2.0) < np_type(1.0) + epsilon}")
print(
f"1 + eps/2 < 1 + eps + eps / 2: {np_type(1.0) + epsilon / np_type(2.0) < np_type(1.0) + epsilon + epsilon / np_type(2.0)}" # noqa: E501
)
print(
f"1 + eps < 1 + eps + eps / 2: {np_type(1.0) + epsilon < np_type(1.0) + epsilon + epsilon / np_type(2.0)}" # noqa: E501
)