Coverage for src/algolib/numerics/constants.py: 100%

1# src/algolib/numerics/constants.py 

2""" 

3Centralized numerical constants for algolib. 

4 

5- No imports from `math` or third-party libs: *pure Python floats only*. 

6- Values are given either as exact decimal, or with comments showing hex-float 

7 for traceability. 

8- Includes Cody–Waite style splits for stable range-reduction in exp/sin/cos. 

9""" 

10 

11from algolib.exceptions import NumericOverflowError 

12 

13# ----------------------------- 

14# IEEE-754 double "machine" info 

15# ----------------------------- 

16# Python float is IEEE-754 binary64 on CPython. 

17DBL_MANT_DIG = 53 # significand bits 

18DBL_EPS = 2.0 ** -52 # ~= 2.220446049250313e-16 (unit roundoff) 

19DBL_MIN_EXP = -1022 # min exponent for normal numbers 

20DBL_MAX_EXP = 1024 # max exponent (exclusive bound for 2^e) 

21DBL_MIN = 2.0 ** -1022 # ~= 2.2250738585072014e-308 (min normal) 

22DBL_DENORM_MIN = 2.0 ** -1074 # ~= 5e-324 (min subnormal) 

23DBL_MAX = (2.0 - 2.0 ** -52) * (2.0 ** 1023) # ~= 1.7976931348623157e+308 

24 

25# -------------------------------- 

26# Fundamental mathematical constants 

27# -------------------------------- 

28PI = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862090 

29PI_2 = 1.5707963267948966192313216916397514420985846996875529104874722961539082031431045 # π/2 

30PI_4 = 0.78539816339744830961566084581987572104929234984377645524373614807695410157155225 # π/4 

31TAU = 6.2831853071795864769252867665590057683943387987502116419498891846156328125724180 # 2π 

32INV_PI = 3.1830988618379067153776752674502872406891929148091289749533468811779359526845307e-1 # 1/π 

33INV_PI_2 = 6.3661977236758134307553505349005744813783858296182579499066937623558719053690614e-1 # 2/π 

34E = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535476 

35LN2 = 6.9314718055994530941723212145817656807550013436025525412068000949339362196969472e-1 

36INV_LN2 = 1.4426950408889634073599246810018921374266459541529859341354494069311092191811851 # 1/ln 2 

37LOG2E = INV_LN2 # alias 

38LN10 = 2.3025850929940456840179914546843642076011014886287729760333279009675726096773525 

39INV_LN10 = 4.3429448190325182765112891891660508229439700580366656611445378316586464920887077e-1 

40SQRT2 = 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070 

41INV_SQRT2 = 7.0710678118654752440084436210484903928483593768847403658833986899536623923105352e-1 

42 

43# ---------------------------------------------------- 

44# Cody–Waite splits for stable range reduction / exp() 

45# ---------------------------------------------------- 

46# ln(2) = LN2_HI + LN2_LO (sum exactly equals LN2 in double) 

47LN2_HI = 0.6931471803691238 

48LN2_LO = 1.9082149292705877e-10 

49 

50# π/2 triple split for high-accuracy range reduction 

51# Values compatible with fdlibm-style splits 

52PI2_HI = 1.5707963267948966 

53PI2_MID = 1.9231321691639753e-17 # 0x3c91a62633145c07 

54PI2_LO = -1.5579014153003124e-33 # 0x3b2a19c036e2eb4f 

55 

56# 2/π double split for accurate k = round(x*(2/π)) 

57# High part truncated to ~28 bits precision; low part is the residual. 

58INV_PI_2_HI = 0.6366197723675814 # high part of 2/π 

59INV_PI_2_LO = -5.692446494650995e-17 

60 

61# Optional: π split (rarely needed if you reduce by π/2) 

62PI_HI = 3.141592653589793 

63PI_LO = 1.2246467991473532e-16 

64 

65# ----------------------------------- 

66# Default, *general* numerical tolerances 

67# ----------------------------------- 

68# These are library-wide suggestions; tests may override. 

69REL_EPS_DEFAULT = 2e-12 

70ABS_EPS_DEFAULT = 1e-12 

71 

72# ------------------------------------------------ 

73# Small helpers (no math.*; keep them dependency-free) 

74# ------------------------------------------------ 

75def pow2_int(k: int) -> float: 

76 """Compute 2**k using only multiplies (supports negative k) with IEEE754-style bounds. 

77 

78 Behavior: 

79 - k > 1023 -> raise NumericOverflowError (subclass of OverflowError) 

80 - k < -1074 -> 0.0 (underflow to zero) 

81 - -1074 <= k <= -1023 -> exact subnormal via halving from DBL_MIN 

82 - otherwise -> exponentiation-by-squaring (<= 1023 以及 >= -1022) 

83 """ 

84 if k == 0: 

85 return 1.0 

86 if k > 1023: 

87 # mirror builtin overflow, 但用库内异常（是 OverflowError 的子类） 

88 raise NumericOverflowError("Result too large") 

89 if k < -1074: 

90 # underflow to zero 

91 return 0.0 

92 

93 if k >= DBL_MIN_EXP: # 正规区间（包含 -1022） 

94 if k > 0: 

95 y = 1.0 

96 base = 2.0 

97 n = k 

98 while n: 

99 if (n & 1) != 0: 

100 y *= base 

101 base *= base 

102 n >>= 1 

103 return y 

104 else: 

105 # -1022 <= k <= -1：用 0.5 的平方幂，稳定且不会一路掉成 0 

106 y = 1.0 

107 base = 0.5 # 2**-1 

108 n = -k 

109 while n: 

110 if (n & 1) != 0: 

111 y *= base 

112 base *= base 

113 n >>= 1 

114 return y 

115 

116 # 次正规：-1074 <= k <= -1023 

117 # 思路：从最小正规 DBL_MIN = 2**-1022 开始，做 (-(k+1022)) 次“乘 0.5”。 

118 # 这里最多 52 次，能得到精确的 subnormal。 

119 steps = -(k + 1022) # 介于 1..52 

120 y = DBL_MIN 

121 for _ in range(steps): 

122 y *= 0.5 

123 # k = -1074 时 y == DBL_DENORM_MIN（5e-324） 

124 return y 

125 

126def copysign1(x: float, y: float) -> float: 

127 r"""Return :math:\\abs{x} with the sign of y. No math module, handles ±0.0 and NaN.""" 

128 # y 是 NaN：按约定返回 x 本身（测试也是这么期望的） 

129 if y != y: # NaN 自不等 

130 return x 

131 

132 # 取 |x| 

133 ax = x if x >= 0.0 else -x 

134 

135 # y == 0.0 时，用十六进制表示判断符号位 

136 if y == 0.0: 

137 neg = float.hex(y).startswith("-") 

138 return -ax if neg else ax 

139 

140 # 其余情况：比较即可 

141 return ax if y > 0.0 else -ax