Coverage for moptipy/utils/math.py: 85%

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.

90 >>> try:

91 ... try_int(-inf)

92 ... except ValueError as ve:

93 ... print(ve)

94 Value must be finite, but is -inf.

96 >>> try:

97 ... try_int(nan)

98 ... except ValueError as ve:

99 ... print(ve)

100 Value must be finite, but is nan.

101

102 >>> try:

103 ... try_int("blab") # noqa # type: off

104 ... except TypeError as te:

105 ... print(te)

106 val should be an instance of any in {float, int} but is str, namely \

107'blab'.

108

109 >>> type(try_int(9007199254740992.0))

110 <class 'int'>

111 >>> try_int(9007199254740992.0)

112 9007199254740992

113 >>> too_big = nextafter(9007199254740992.0, inf)

114 >>> print(too_big)

115 9007199254740994.0

116 >>> type(try_int(too_big))

117 <class 'float'>

118 >>> type(try_int(-9007199254740992.0))

119 <class 'int'>

120 >>> try_int(-9007199254740992.0)

121 -9007199254740992

122 >>> type(try_int(-too_big))

123 <class 'float'>

124 """

125 if isinstance(val, int):

126 return val

127 if isinstance(val, float):

128 if not isfinite(val):

129 raise ValueError(f"Value must be finite, but is {val}.")

130 if _DBL_INT_LIMIT_N <= val <= DBL_INT_LIMIT_P:

131 a = int(val)

132 if a == val:

133 return a

134 return val

135 raise type_error(val, "val", (int, float))

136

137

138def try_int_div(a: int, b: int) -> int | float:

139 """

140 Try to divide two integers at best precision.

141

142 Floating point divisions can incur some loss of precision. We try

143 to avoid this here as much as possible. First, we check if `a` is

144 divisible by `b` without any fractional part. If this is true, then

145 we can do a pure integer division without loss of any precision.

146

147 Otherwise, we would do a floating point division. This will lead the

148 values be converted to floating point numbers, which can incur

149 a loss of precision. In the past, we tried to divide both numbers by

150 their `gcd` to make them smaller in order to avoid a loss of precision.

151 But it seems that Python is already doing something like this internally

152 when performing the `a / b` division anyway, so we don't do this anymore.

153 However, it is still attempted to convert the result back to an integer.

154

155 :param a: the first integer

156 :param b: the second integer

157 :return: a/b, either as `int` or as `float` but always a finite value

158 :raises OverflowError: if the division must be performed using floating

159 point arithmetic, but the result would be too large to fit into a

160 `float`

161 :raises ZeroDivisionError: if `b==0`

162

163 >>> print(try_int_div(10, 2))

164 5

165 >>> print(try_int_div(10, -2))

166 -5

167 >>> print(try_int_div(-10, 2))

168 -5

169 >>> print(try_int_div(-10, -2))

170 5

171 >>> print(type(try_int_div(10, 2)))

172 <class 'int'>

173 >>> print(try_int_div(10, 3))

174 3.3333333333333335

175 >>> print(try_int_div(-10, 3))

176 -3.3333333333333335

177 >>> print(try_int_div(10, -3))

178 -3.3333333333333335

179 >>> print(try_int_div(-10, -3))

180 3.3333333333333335

181 >>> print(type(try_int_div(10, 3)))

182 <class 'float'>

183 >>> print(try_int_div(9007199254740992, 1))

184 9007199254740992

185 >>> print(try_int_div(2109792310235001520128, 234234))

186 9007199254740992

187 >>> print(try_int_div(2109792310235001520128, 234235))

188 9007160801054503

189 >>> print(try_int_div(2109792310235001520128, 234233))

190 9007237708755818.0

191 >>> large = 123456789012345678901234567890123456789012345678901234567\

19289012345678901234567890123456789012345678901234567890123456789012345678901234\

19356789012345678901234567890123456789012345678901234567890123456789012345678901\

19423456789012345678901234567890123456789012345678901234567890123456789012345678\

19590123456789012345678901234567890123456789012345678901234567890123456789012345\

19667890123456789012345678901234567890123456789012345678901234567890123456789012\

1973456789012345678901234567890123456789012345678901234567890123456789012345678\

19890123456789012345678901234567890123456789012345678901234567890123456789012345\

199678901234567890123456789012345678901234567890

200 >>> try:

201 ... large / 1

202 ... except OverflowError as oe:

203 ... print(oe)

204 integer division result too large for a float

205 >>> try_int_div(large, 1)

206 123456789012345678901234567890123456789012345678901234567\

20789012345678901234567890123456789012345678901234567890123456789012345678901234\

20856789012345678901234567890123456789012345678901234567890123456789012345678901\

20923456789012345678901234567890123456789012345678901234567890123456789012345678\

21090123456789012345678901234567890123456789012345678901234567890123456789012345\

21167890123456789012345678901234567890123456789012345678901234567890123456789012\

2123456789012345678901234567890123456789012345678901234567890123456789012345678\

21390123456789012345678901234567890123456789012345678901234567890123456789012345\

214678901234567890123456789012345678901234567890

215 >>> try_int_div(large * 7, 1 * 7)

216 123456789012345678901234567890123456789012345678901234567\

21789012345678901234567890123456789012345678901234567890123456789012345678901234\

21856789012345678901234567890123456789012345678901234567890123456789012345678901\

21923456789012345678901234567890123456789012345678901234567890123456789012345678\

22090123456789012345678901234567890123456789012345678901234567890123456789012345\

22167890123456789012345678901234567890123456789012345678901234567890123456789012\

2223456789012345678901234567890123456789012345678901234567890123456789012345678\

22390123456789012345678901234567890123456789012345678901234567890123456789012345\

224678901234567890123456789012345678901234567890

225 >>> try:

226 ... try_int_div(large, 7)

227 ... except OverflowError as oe:

228 ... print(oe)

229 integer division result too large for a float

230 >>> try:

231 ... try_int_div(0, 0)

232 ... except ZeroDivisionError as zde:

233 ... print(zde)

234 integer division or modulo by zero

235 >>> try:

236 ... try_int_div(1, 0)

237 ... except ZeroDivisionError as zde:

238 ... print(zde)

239 integer division or modulo by zero

240 >>> try:

241 ... try_int_div(-1, 0)

242 ... except ZeroDivisionError as zde:

243 ... print(zde)

244 integer division or modulo by zero

245 """

246 # First try pure integer division, in case `a` is divisible by `b`.

247 int_div_test: Final[int] = a // b

248 if (int_div_test * b) == a:

249 return int_div_test

250 # If that does not work, use the floating point division.

251 return __try_int(a / b)

252

253

254def try_float_div(a: int | float, b: int | float) -> int | float:

255 """

256 Try to divide two numbers at best precision.

257

258 First, we will check if we can convert the numbers to integers

259 without loss of precision via :func:`try_int`. If yes, then

260 we go for the maximum-precision integer division via :func:`try_int_div`.

261 If no, then we do the normal floating point division and try to convert

262 the result to an integer if that can be done without loss of precision.

263

264 :param a: the first number

265 :param b: the second number

266 :return: `a/b`, but always finite

267

268 :raises ValueError: if either one of the arguments or the final result

269 would not be finite

270

271 >>> try_float_div(1e180, 1e60)

272 1.0000000000000001e+120

273 >>> try_float_div(1e60, 1e-60)

274 1e+120

275 >>> try_float_div(1e14, 1e-1)

276 1000000000000000

277 >>> try_float_div(1e14, -1e-1)

278 -1000000000000000

279 >>> try_float_div(-1e14, 1e-1)

280 -1000000000000000

281 >>> try_float_div(-1e14, -1e-1)

282 1000000000000000

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-1)

288 -1e+16

289 >>> try_float_div(-1e15, -1e-1)

290 1e+16

291 >>> try_float_div(1e15, 1e-15)

292 9.999999999999999e+29

293

294 >>> print(type(try_float_div(10, 2)))

295 <class 'int'>

296 >>> print(type(try_float_div(10, 3)))

297 <class 'float'>

298 >>> print(type(try_float_div(10, 0.5)))

299 <class 'int'>

300

301 >>> from math import inf, nan

302 >>> try:

303 ... try_float_div(1.0, 0)

304 ... except ZeroDivisionError as zde:

305 ... print(zde)

306 integer division or modulo by zero

307

308 >>> try:

309 ... try_float_div(1.0, -0.0)

310 ... except ZeroDivisionError as zde:

311 ... print(zde)

312 integer division or modulo by zero

313

314 >>> try:

315 ... try_float_div(inf, 0)

316 ... except ValueError as ve:

317 ... print(ve)

318 Value must be finite, but is inf.

319

320 >>> try:

321 ... try_float_div(-inf, 0)

322 ... except ValueError as ve:

323 ... print(ve)

324 Value must be finite, but is -inf.

325

326 >>> try:

327 ... try_float_div(nan, 0)

328 ... except ValueError as ve:

329 ... print(ve)

330 Value must be finite, but is nan.

331

332 >>> try:

333 ... try_float_div(1, inf)

334 ... except ValueError as ve:

335 ... print(ve)

336 Value must be finite, but is inf.

337

338 >>> try:

339 ... try_float_div(1, -inf)

340 ... except ValueError as ve:

341 ... print(ve)

342 Value must be finite, but is -inf.

343

344 >>> try:

345 ... try_float_div(1, nan)

346 ... except ValueError as ve:

347 ... print(ve)

348 Value must be finite, but is nan.

349

350 >>> try:

351 ... try_float_div(1e300, 1e-60)

352 ... except ValueError as ve:

353 ... print(ve)

354 Result must be finite, but is 1e+300/1e-60=inf.

355 """

356 ia: Final[int | float] = try_int(a)

357 ib: Final[int | float] = try_int(b)

358 if isinstance(ia, int) and isinstance(ib, int):

359 return try_int_div(ia, ib)

360 val: Final[float] = ia / ib

361 if not isfinite(val):

362 raise ValueError(f"Result must be finite, but is {a}/{b}={val}.")

363 return __try_int(val)

364

365

366def __bin_search_root(value: int, power: int) -> int:

367 """

368 Search a truncated integer root via binary search.

369

370 :param value: the integer value to compute the root of

371 :param power: the integer power of the root

372 :returns: the integer root, truncated if necessary

373 """

374 int_sqrt: Final[int] = isqrt(value)

375 if power == 2:

376 return int_sqrt

377

378 # compute the maximum base root

379 root_min: int = 1

380 root_max: int = (int_sqrt - 1) if value > 3 else 1

381 while root_max >= root_min:

382 root_mid = isqrt(root_min * root_max)

383 power_mid = root_mid ** power

384 if power_mid > value:

385 root_max = root_mid - 1

386 elif power_mid < value:

387 root_min = root_mid + 1

388 else:

389 return root_mid

390 return root_max

391

392

393def try_int_root(value: int, power: int,

394 none_on_overflow: bool = True) -> int | float | None:

395 """

396 Compute `value**(1/power)` where `value` and `power` are both integers.

397

398 :param value: the integer value to compute the root of

399 :param power: the integer power of the root

400 :param none_on_overflow: `True` if `None` should be returned if we

401 encounter an :class:`OverflowError`, `False` if we should simply

402 re-raise it

403 :returns: the root, or `None` if we encounter an overflow

404

405 >>> try_int_root(100, 2)

406 10

407 >>> try_int_root(100, 3) - (100 ** (1/3))

408 0.0

409 >>> try_int_root(123100, 23) - (123100 ** (1/23))

410 0.0

411 >>> 123100**1000 - try_int_root(123100**1000, 1000)**1000

412 0

413 >>> 11**290 - try_int_root(11**290, 290)**290

414 0

415 >>> (abs(1.797E308 - try_int_root(int(1.797E308), 10) ** 10)

416 ... < abs(1.797E308 - (int(1.797E308) ** 0.1) ** 10))

417 True

418 """

419 if not isinstance(value, int):

420 raise type_error(value, "value", int)

421 if not isinstance(power, int):

422 raise type_error(power, "power", int)

423 if power <= 0:

424 raise ValueError(f"power must be positive but is {power}.")

425 if value <= 1:

426 if value < 0:

427 raise ValueError(f"value must not be negative, but is {value}.")

428 return value

429 if power == 1:

430 return value

431

432 # first, approximate root with crude binary search

433 int_root: int = __bin_search_root(value, power)

434 root_power: int = int_root ** power

435 if root_power >= value:

436 if root_power == value:

437 return int_root # ok, we are done

438 raise ValueError(f"{int_root}**{power}={root_power} > {value}?")

439

440 # Now try to remove integer factors from the value and the root,

441 # to achieve a more accurate division of the rest.

442 # This is probably extremely inefficient.

443 root_base: int = 1

444 i: int = 2

445 end: int = min(10_000, int_root)

446 while i < end:

447 div: int = i ** power

448 if (value % div) == 0:

449 root_base *= i

450 value //= div

451 int_root = (int_root + i - 1) // i

452 end = min(end, int_root)

453 else:

454 i += 1

455

456 # value is now reduced by any integer factor that we can pull out

457 # of the root. These factors are aggregated in root_base.

458 # If root_base != 1, we need to re-compute the root of the remainder

459 # inside value. Whatever root we now compute of value, we must multiply

460 # it with root_base.

461 if root_base != 1:

462 int_root = __bin_search_root(value, power)

463

464 # check again

465 root_power = int_root ** power

466 if root_power >= value:

467 if root_power == value:

468 return root_base * int_root # ok, we are done

469 raise ValueError(f"{int_root}**{power}={root_power} > {value}?")

470

471 # from now on, root may be either and int or (more likely) a float.

472 root: int | float = int_root

473 try:

474 rest: int | float = try_int_div(value, root_power)

475 rest_root: int | float = __try_int(rest ** (1.0 / power))

476 root = __try_int(root * rest_root)

477 except OverflowError: # as ofe:

478 if none_on_overflow:

479 return None

480 raise # raises ofe again

481

482 # OK, we got an approximate root of what remains of value.

483 # Let's see if we can refine it.

484 with contextlib.suppress(OverflowError):

485 diff = abs((root ** power) - value)

486 root2: int | float = __try_int(exp(log(value) / root))

487 diff2: int | float = abs((root2 ** power) - value)

488 if diff2 < diff:

489 diff = diff2

490 root = root2

491 root2 = __try_int(2 ** (log2(value) / root))

492 diff2 = abs((root2 ** power) - value)

493 if diff2 < diff:

494 diff = diff2

495 root = root2

496

497 rdn = root

498 rup = root

499 while True:

500 rdn = nextafter(rdn, -inf)

501 apd = abs((rdn ** power) - value)

502 if apd > diff:

503 break

504 diff = apd

505 root = rdn

506 while True:

507 rup = nextafter(rup, inf)

508 apd = abs((rup ** power) - value)

509 if apd > diff:

510 break

511 diff = apd

512 root = rup

513

514 return root_base * __try_int(root)

Coverage for moptipy / utils / math.py: 85%

128 statements