@@ -21,37 +21,14 @@ int64_t gcd64(int64_t a, int64_t b) {
2121 return a;
2222}
2323
24- __int128 gcd128 (__int128 a, __int128 b) {
25- if (a < 0 ) a = -a;
26- if (b < 0 ) b = -b;
27- while (b) {
28- __int128 t = a % b;
29- a = b;
30- b = t;
31- }
32- return a;
24+ // Overflow-checked 64-bit arithmetic. The target (32-bit Cortex-M7) has no
25+ // __int128, so we rely on the compiler builtins and demote to double whenever
26+ // an exact result would overflow int64.
27+ inline bool mulOverflow (int64_t a, int64_t b, int64_t * out) {
28+ return __builtin_mul_overflow (a, b, out);
3329}
34-
35- bool fitsInt64 (__int128 x) {
36- return x <= (__int128)INT64_MAX && x >= (__int128)INT64_MIN ;
37- }
38-
39- // Build an exact value from a 128-bit fraction, demoting to APPROX on overflow.
40- Value fromFraction128 (__int128 num, __int128 den) {
41- if (den == 0 ) return Value::undefined ();
42- if (den < 0 ) {
43- num = -num;
44- den = -den;
45- }
46- __int128 g = gcd128 (num, den);
47- if (g != 0 ) {
48- num /= g;
49- den /= g;
50- }
51- if (fitsInt64 (num) && fitsInt64 (den)) {
52- return Value::rational ((int64_t )num, (int64_t )den);
53- }
54- return Value::real ((double )num / (double )den);
30+ inline bool addOverflow (int64_t a, int64_t b, int64_t * out) {
31+ return __builtin_add_overflow (a, b, out);
5532}
5633
5734int64_t isqrt64 (int64_t n) {
@@ -100,9 +77,13 @@ Value Value::inverse() const {
10077Value Value::add (const Value& a, const Value& b) {
10178 if (a.m_undef || b.m_undef ) return undefined ();
10279 if (a.isExact () && b.isExact ()) {
103- return fromFraction128 (
104- (__int128)a.m_num * b.m_den + (__int128)b.m_num * a.m_den ,
105- (__int128)a.m_den * b.m_den );
80+ int64_t left, right, num, den;
81+ if (!mulOverflow (a.m_num , b.m_den , &left) &&
82+ !mulOverflow (b.m_num , a.m_den , &right) &&
83+ !addOverflow (left, right, &num) &&
84+ !mulOverflow (a.m_den , b.m_den , &den)) {
85+ return rational (num, den);
86+ }
10687 }
10788 return real (a.toDouble () + b.toDouble ());
10889}
@@ -114,8 +95,11 @@ Value Value::sub(const Value& a, const Value& b) {
11495Value Value::mul (const Value& a, const Value& b) {
11596 if (a.m_undef || b.m_undef ) return undefined ();
11697 if (a.isExact () && b.isExact ()) {
117- return fromFraction128 ((__int128)a.m_num * b.m_num ,
118- (__int128)a.m_den * b.m_den );
98+ int64_t num, den;
99+ if (!mulOverflow (a.m_num , b.m_num , &num) &&
100+ !mulOverflow (a.m_den , b.m_den , &den)) {
101+ return rational (num, den);
102+ }
119103 }
120104 return real (a.toDouble () * b.toDouble ());
121105}
@@ -124,8 +108,11 @@ Value Value::div(const Value& a, const Value& b) {
124108 if (a.m_undef || b.m_undef ) return undefined ();
125109 if (a.isExact () && b.isExact ()) {
126110 if (b.m_num == 0 ) return undefined ();
127- return fromFraction128 ((__int128)a.m_num * b.m_den ,
128- (__int128)a.m_den * b.m_num );
111+ int64_t num, den;
112+ if (!mulOverflow (a.m_num , b.m_den , &num) &&
113+ !mulOverflow (a.m_den , b.m_num , &den)) {
114+ return rational (num, den);
115+ }
129116 }
130117 double d = b.toDouble ();
131118 if (d == 0.0 ) return undefined ();
@@ -134,22 +121,23 @@ Value Value::div(const Value& a, const Value& b) {
134121
135122Value Value::pow (const Value& a, const Value& b) {
136123 if (a.m_undef || b.m_undef ) return undefined ();
137- // Exact when the exponent is an integer that fits a reasonable range.
124+ // Exact when the exponent is an integer in a reasonable range.
138125 if (a.isExact () && b.isExact () && b.m_den == 1 && b.m_num <= 63 &&
139126 b.m_num >= -63 ) {
140127 int64_t e = b.m_num ;
141128 bool negativeExponent = e < 0 ;
142129 if (negativeExponent) e = -e;
143- __int128 num = 1 , den = 1 ;
144- for (int64_t i = 0 ; i < e; i++) {
145- num *= a.m_num ;
146- den *= a.m_den ;
147- if (!fitsInt64 (num) || !fitsInt64 (den)) {
148- return real (::pow (a.toDouble (), b.toDouble ()));
149- }
130+ int64_t num = 1 , den = 1 ;
131+ bool overflow = false ;
132+ for (int64_t i = 0 ; i < e && !overflow; i++) {
133+ overflow = mulOverflow (num, a.m_num , &num) ||
134+ mulOverflow (den, a.m_den , &den);
135+ }
136+ if (!overflow) {
137+ Value result = rational (num, den);
138+ return negativeExponent ? result.inverse () : result;
150139 }
151- Value result = fromFraction128 (num, den);
152- return negativeExponent ? result.inverse () : result;
140+ return real (::pow (a.toDouble (), b.toDouble ()));
153141 }
154142 double base = a.toDouble ();
155143 double exp = b.toDouble ();
0 commit comments