Coverage for moptipy / utils / math.py: 85%
128 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-11-24 08:49 +0000
1"""Simple routines for handling numbers and numerical stuff."""
3import contextlib
4from math import exp, inf, isfinite, isqrt, log, log2, nextafter
5from typing import Final
7from pycommons.types import type_error
9#: The positive limit for doubles that can be represented exactly as ints.
10DBL_INT_LIMIT_P: Final[float] = 9007199254740992.0 # = 1 << 53
11#: The negative limit for doubles that can be represented exactly as ints.
12_DBL_INT_LIMIT_N: Final[float] = -DBL_INT_LIMIT_P
15def __try_int(val: float) -> int | float:
16 """
17 Convert a float to an int without any fancy checks.
19 :param val: the flot
20 :returns: the float or int
22 >>> from math import inf, nan, nextafter
23 >>> type(__try_int(0.0))
24 <class 'int'>
25 >>> type(__try_int(0.5))
26 <class 'float'>
27 >>> type(__try_int(inf))
28 <class 'float'>
29 >>> type(__try_int(-inf))
30 <class 'float'>
31 >>> type(__try_int(nan))
32 <class 'float'>
33 >>> 1 << 53
34 9007199254740992
35 >>> type(__try_int(9007199254740992.0))
36 <class 'int'>
37 >>> __try_int(9007199254740992.0)
38 9007199254740992
39 >>> too_big = nextafter(9007199254740992.0, inf)
40 >>> print(too_big)
41 9007199254740994.0
42 >>> type(__try_int(too_big))
43 <class 'float'>
44 >>> type(__try_int(-9007199254740992.0))
45 <class 'int'>
46 >>> __try_int(-9007199254740992.0)
47 -9007199254740992
48 >>> type(__try_int(-too_big))
49 <class 'float'>
50 """
51 if _DBL_INT_LIMIT_N <= val <= DBL_INT_LIMIT_P:
52 a = int(val)
53 if a == val:
54 return a
55 return val
58def try_int(val: int | float) -> int | float:
59 """
60 Attempt to convert a float to an integer.
62 This method will convert a floating point number to an integer if the
63 floating point number was representing an exact integer. This is the
64 case if it has a) no fractional part and b) is in the range
65 `-9007199254740992...9007199254740992`, i.e., the range where `+1` and
66 `-1` work without loss of precision.
68 :param val: the input value, which must either be `int` or `float`
69 :return: an `int` if `val` can be represented as `int` without loss of
70 precision, `val` otherwise
71 :raises TypeError: if `val` is neither an instance of `int` nor of `float`
72 :raises ValueError: if `val` is a `float`, but not finite
74 >>> print(type(try_int(10.5)))
75 <class 'float'>
76 >>> print(type(try_int(10)))
77 <class 'int'>
79 >>> from math import inf, nan, nextafter
80 >>> type(try_int(0.0))
81 <class 'int'>
82 >>> type(try_int(0.5))
83 <class 'float'>
84 >>> try:
85 ... try_int(inf)
86 ... except ValueError as ve:
87 ... print(ve)
88 Value must be finite, but is inf.
89 >>> try:
90 ... try_int(-inf)
91 ... except ValueError as ve:
92 ... print(ve)
93 Value must be finite, but is -inf.
94 >>> try:
95 ... try_int(nan)
96 ... except ValueError as ve:
97 ... print(ve)
98 Value must be finite, but is nan.
99 >>> try:
100 ... try_int("blab") # noqa # type: off
101 ... except TypeError as te:
102 ... print(te)
103 val should be an instance of any in {float, int} but is str, namely \
104'blab'.
105 >>> type(try_int(9007199254740992.0))
106 <class 'int'>
107 >>> try_int(9007199254740992.0)
108 9007199254740992
109 >>> too_big = nextafter(9007199254740992.0, inf)
110 >>> print(too_big)
111 9007199254740994.0
112 >>> type(try_int(too_big))
113 <class 'float'>
114 >>> type(try_int(-9007199254740992.0))
115 <class 'int'>
116 >>> try_int(-9007199254740992.0)
117 -9007199254740992
118 >>> type(try_int(-too_big))
119 <class 'float'>
120 """
121 if isinstance(val, int):
122 return val
123 if isinstance(val, float):
124 if not isfinite(val):
125 raise ValueError(f"Value must be finite, but is {val}.")
126 if _DBL_INT_LIMIT_N <= val <= DBL_INT_LIMIT_P:
127 a = int(val)
128 if a == val:
129 return a
130 return val
131 raise type_error(val, "val", (int, float))
134def try_int_div(a: int, b: int) -> int | float:
135 """
136 Try to divide two integers at best precision.
138 Floating point divisions can incur some loss of precision. We try
139 to avoid this here as much as possible. First, we check if `a` is
140 divisible by `b` without any fractional part. If this is true, then
141 we can do a pure integer division without loss of any precision.
143 Otherwise, we would do a floating point division. This will lead the
144 values be converted to floating point numbers, which can incur
145 a loss of precision. In the past, we tried to divide both numbers by
146 their `gcd` to make them smaller in order to avoid a loss of precision.
147 But it seems that Python is already doing something like this internally
148 when performing the `a / b` division anyway, so we don't do this anymore.
149 However, it is still attempted to convert the result back to an integer.
151 :param a: the first integer
152 :param b: the second integer
153 :return: a/b, either as `int` or as `float` but always a finite value
154 :raises OverflowError: if the division must be performed using floating
155 point arithmetic, but the result would be too large to fit into a
156 `float`
157 :raises ZeroDivisionError: if `b==0`
159 >>> print(try_int_div(10, 2))
160 5
161 >>> print(try_int_div(10, -2))
162 -5
163 >>> print(try_int_div(-10, 2))
164 -5
165 >>> print(try_int_div(-10, -2))
166 5
167 >>> print(type(try_int_div(10, 2)))
168 <class 'int'>
169 >>> print(try_int_div(10, 3))
170 3.3333333333333335
171 >>> print(try_int_div(-10, 3))
172 -3.3333333333333335
173 >>> print(try_int_div(10, -3))
174 -3.3333333333333335
175 >>> print(try_int_div(-10, -3))
176 3.3333333333333335
177 >>> print(type(try_int_div(10, 3)))
178 <class 'float'>
179 >>> print(try_int_div(9007199254740992, 1))
180 9007199254740992
181 >>> print(try_int_div(2109792310235001520128, 234234))
182 9007199254740992
183 >>> print(try_int_div(2109792310235001520128, 234235))
184 9007160801054503
185 >>> print(try_int_div(2109792310235001520128, 234233))
186 9007237708755818.0
187 >>> large = 123456789012345678901234567890123456789012345678901234567\
18889012345678901234567890123456789012345678901234567890123456789012345678901234\
18956789012345678901234567890123456789012345678901234567890123456789012345678901\
19023456789012345678901234567890123456789012345678901234567890123456789012345678\
19190123456789012345678901234567890123456789012345678901234567890123456789012345\
19267890123456789012345678901234567890123456789012345678901234567890123456789012\
1933456789012345678901234567890123456789012345678901234567890123456789012345678\
19490123456789012345678901234567890123456789012345678901234567890123456789012345\
195678901234567890123456789012345678901234567890
196 >>> try:
197 ... large / 1
198 ... except OverflowError as oe:
199 ... print(oe)
200 integer division result too large for a float
201 >>> try_int_div(large, 1)
202 123456789012345678901234567890123456789012345678901234567\
20389012345678901234567890123456789012345678901234567890123456789012345678901234\
20456789012345678901234567890123456789012345678901234567890123456789012345678901\
20523456789012345678901234567890123456789012345678901234567890123456789012345678\
20690123456789012345678901234567890123456789012345678901234567890123456789012345\
20767890123456789012345678901234567890123456789012345678901234567890123456789012\
2083456789012345678901234567890123456789012345678901234567890123456789012345678\
20990123456789012345678901234567890123456789012345678901234567890123456789012345\
210678901234567890123456789012345678901234567890
211 >>> try_int_div(large * 7, 1 * 7)
212 123456789012345678901234567890123456789012345678901234567\
21389012345678901234567890123456789012345678901234567890123456789012345678901234\
21456789012345678901234567890123456789012345678901234567890123456789012345678901\
21523456789012345678901234567890123456789012345678901234567890123456789012345678\
21690123456789012345678901234567890123456789012345678901234567890123456789012345\
21767890123456789012345678901234567890123456789012345678901234567890123456789012\
2183456789012345678901234567890123456789012345678901234567890123456789012345678\
21990123456789012345678901234567890123456789012345678901234567890123456789012345\
220678901234567890123456789012345678901234567890
221 >>> try:
222 ... try_int_div(large, 7)
223 ... except OverflowError as oe:
224 ... print(oe)
225 integer division result too large for a float
226 >>> try:
227 ... try_int_div(0, 0)
228 ... except ZeroDivisionError as zde:
229 ... print(zde)
230 integer division or modulo by zero
231 >>> try:
232 ... try_int_div(1, 0)
233 ... except ZeroDivisionError as zde:
234 ... print(zde)
235 integer division or modulo by zero
236 >>> try:
237 ... try_int_div(-1, 0)
238 ... except ZeroDivisionError as zde:
239 ... print(zde)
240 integer division or modulo by zero
241 """
242 # First try pure integer division, in case `a` is divisible by `b`.
243 int_div_test: Final[int] = a // b
244 if (int_div_test * b) == a:
245 return int_div_test
246 # If that does not work, use the floating point division.
247 return __try_int(a / b)
250def try_float_div(a: int | float, b: int | float) -> int | float:
251 """
252 Try to divide two numbers at best precision.
254 First, we will check if we can convert the numbers to integers
255 without loss of precision via :func:`try_int`. If yes, then
256 we go for the maximum-precision integer division via :func:`try_int_div`.
257 If no, then we do the normal floating point division and try to convert
258 the result to an integer if that can be done without loss of precision.
260 :param a: the first number
261 :param b: the second number
262 :return: `a/b`, but always finite
264 :raises ValueError: if either one of the arguments or the final result
265 would not be finite
267 >>> try_float_div(1e180, 1e60)
268 1.0000000000000001e+120
269 >>> try_float_div(1e60, 1e-60)
270 1e+120
271 >>> try_float_div(1e14, 1e-1)
272 1000000000000000
273 >>> try_float_div(1e14, -1e-1)
274 -1000000000000000
275 >>> try_float_div(-1e14, 1e-1)
276 -1000000000000000
277 >>> try_float_div(-1e14, -1e-1)
278 1000000000000000
279 >>> try_float_div(1e15, 1e-1)
280 1e+16
281 >>> try_float_div(1e15, -1e-1)
282 -1e+16
283 >>> try_float_div(-1e15, 1e-1)
284 -1e+16
285 >>> try_float_div(-1e15, -1e-1)
286 1e+16
287 >>> try_float_div(1e15, 1e-15)
288 9.999999999999999e+29
290 >>> print(type(try_float_div(10, 2)))
291 <class 'int'>
292 >>> print(type(try_float_div(10, 3)))
293 <class 'float'>
294 >>> print(type(try_float_div(10, 0.5)))
295 <class 'int'>
297 >>> from math import inf, nan
298 >>> try:
299 ... try_float_div(1.0, 0)
300 ... except ZeroDivisionError as zde:
301 ... print(zde)
302 integer division or modulo by zero
304 >>> try:
305 ... try_float_div(1.0, -0.0)
306 ... except ZeroDivisionError as zde:
307 ... print(zde)
308 integer division or modulo by zero
310 >>> try:
311 ... try_float_div(inf, 0)
312 ... except ValueError as ve:
313 ... print(ve)
314 Value must be finite, but is inf.
316 >>> try:
317 ... try_float_div(-inf, 0)
318 ... except ValueError as ve:
319 ... print(ve)
320 Value must be finite, but is -inf.
322 >>> try:
323 ... try_float_div(nan, 0)
324 ... except ValueError as ve:
325 ... print(ve)
326 Value must be finite, but is nan.
328 >>> try:
329 ... try_float_div(1, inf)
330 ... except ValueError as ve:
331 ... print(ve)
332 Value must be finite, but is inf.
334 >>> try:
335 ... try_float_div(1, -inf)
336 ... except ValueError as ve:
337 ... print(ve)
338 Value must be finite, but is -inf.
340 >>> try:
341 ... try_float_div(1, nan)
342 ... except ValueError as ve:
343 ... print(ve)
344 Value must be finite, but is nan.
346 >>> try:
347 ... try_float_div(1e300, 1e-60)
348 ... except ValueError as ve:
349 ... print(ve)
350 Result must be finite, but is 1e+300/1e-60=inf.
351 """
352 ia: Final[int | float] = try_int(a)
353 ib: Final[int | float] = try_int(b)
354 if isinstance(ia, int) and isinstance(ib, int):
355 return try_int_div(ia, ib)
356 val: Final[float] = ia / ib
357 if not isfinite(val):
358 raise ValueError(f"Result must be finite, but is {a}/{b}={val}.")
359 return __try_int(val)
362def __bin_search_root(value: int, power: int) -> int:
363 """
364 Search a truncated integer root via binary search.
366 :param value: the integer value to compute the root of
367 :param power: the integer power of the root
368 :returns: the integer root, truncated if necessary
369 """
370 int_sqrt: Final[int] = isqrt(value)
371 if power == 2:
372 return int_sqrt
374 # compute the maximum base root
375 root_min: int = 1
376 root_max: int = (int_sqrt - 1) if value > 3 else 1
377 while root_max >= root_min:
378 root_mid = isqrt(root_min * root_max)
379 power_mid = root_mid ** power
380 if power_mid > value:
381 root_max = root_mid - 1
382 elif power_mid < value:
383 root_min = root_mid + 1
384 else:
385 return root_mid
386 return root_max
389def try_int_root(value: int, power: int,
390 none_on_overflow: bool = True) -> int | float | None:
391 """
392 Compute `value**(1/power)` where `value` and `power` are both integers.
394 :param value: the integer value to compute the root of
395 :param power: the integer power of the root
396 :param none_on_overflow: `True` if `None` should be returned if we
397 encounter an :class:`OverflowError`, `False` if we should simply
398 re-raise it
399 :returns: the root, or `None` if we encounter an overflow
401 >>> try_int_root(100, 2)
402 10
403 >>> try_int_root(100, 3) - (100 ** (1/3))
404 0.0
405 >>> try_int_root(123100, 23) - (123100 ** (1/23))
406 0.0
407 >>> 123100**1000 - try_int_root(123100**1000, 1000)**1000
408 0
409 >>> 11**290 - try_int_root(11**290, 290)**290
410 0
411 >>> (abs(1.797E308 - try_int_root(int(1.797E308), 10) ** 10)
412 ... < abs(1.797E308 - (int(1.797E308) ** 0.1) ** 10))
413 True
414 """
415 if not isinstance(value, int):
416 raise type_error(value, "value", int)
417 if not isinstance(power, int):
418 raise type_error(power, "power", int)
419 if power <= 0:
420 raise ValueError(f"power must be positive but is {power}.")
421 if value <= 1:
422 if value < 0:
423 raise ValueError(f"value must not be negative, but is {value}.")
424 return value
425 if power == 1:
426 return value
428 # first, approximate root with crude binary search
429 int_root: int = __bin_search_root(value, power)
430 root_power: int = int_root ** power
431 if root_power >= value:
432 if root_power == value:
433 return int_root # ok, we are done
434 raise ValueError(f"{int_root}**{power}={root_power} > {value}?")
436 # Now try to remove integer factors from the value and the root,
437 # to achieve a more accurate division of the rest.
438 # This is probably extremely inefficient.
439 root_base: int = 1
440 i: int = 2
441 end: int = min(10_000, int_root)
442 while i < end:
443 div: int = i ** power
444 if (value % div) == 0:
445 root_base *= i
446 value //= div
447 int_root = (int_root + i - 1) // i
448 end = min(end, int_root)
449 else:
450 i += 1
452 # value is now reduced by any integer factor that we can pull out
453 # of the root. These factors are aggregated in root_base.
454 # If root_base != 1, we need to re-compute the root of the remainder
455 # inside value. Whatever root we now compute of value, we must multiply
456 # it with root_base.
457 if root_base != 1:
458 int_root = __bin_search_root(value, power)
460 # check again
461 root_power = int_root ** power
462 if root_power >= value:
463 if root_power == value:
464 return root_base * int_root # ok, we are done
465 raise ValueError(f"{int_root}**{power}={root_power} > {value}?")
467 # from now on, root may be either and int or (more likely) a float.
468 root: int | float = int_root
469 try:
470 rest: int | float = try_int_div(value, root_power)
471 rest_root: int | float = __try_int(rest ** (1.0 / power))
472 root = __try_int(root * rest_root)
473 except OverflowError: # as ofe:
474 if none_on_overflow:
475 return None
476 raise # raises ofe again
478 # OK, we got an approximate root of what remains of value.
479 # Let's see if we can refine it.
480 with contextlib.suppress(OverflowError):
481 diff = abs((root ** power) - value)
482 root2: int | float = __try_int(exp(log(value) / root))
483 diff2: int | float = abs((root2 ** power) - value)
484 if diff2 < diff:
485 diff = diff2
486 root = root2
487 root2 = __try_int(2 ** (log2(value) / root))
488 diff2 = abs((root2 ** power) - value)
489 if diff2 < diff:
490 diff = diff2
491 root = root2
493 rdn = root
494 rup = root
495 while True:
496 rdn = nextafter(rdn, -inf)
497 apd = abs((rdn ** power) - value)
498 if apd > diff:
499 break
500 diff = apd
501 root = rdn
502 while True:
503 rup = nextafter(rup, inf)
504 apd = abs((rup ** power) - value)
505 if apd > diff:
506 break
507 diff = apd
508 root = rup
510 return root_base * __try_int(root)