Skip to content

Commit e4cd5ee

Browse files
implement pow64_overflow, lose the macro
1 parent 3bf23a6 commit e4cd5ee

2 files changed

Lines changed: 49 additions & 22 deletions

File tree

src/integer64.c

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -366,10 +366,11 @@ SEXP power_integer64_integer64(SEXP e1_, SEXP e2_, SEXP ret_){
366366
long long * e1 = (long long *) REAL(e1_);
367367
long long * e2 = (long long *) REAL(e2_);
368368
long long * ret = (long long *) REAL(ret_);
369-
long long base, exp;
370369
Rboolean naflag = FALSE;
371370
mod_iterate(n1, n2, i1, i2) {
372-
POW64(e1[i1],e2[i2],ret[i],naflag,base,exp)
371+
if (pow64_overflow(e1[i1], e2[i2], &ret[i])) {
372+
naflag = TRUE;
373+
}
373374
}
374375
if (naflag)warning(INTEGER64_OVERFLOW_WARNING);
375376
return ret_;

src/integer64.h

Lines changed: 46 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -118,6 +118,52 @@ static inline bool mul64_overflow(long long a, long long b, long long *res) {
118118
#endif
119119
}
120120

121+
static inline bool pow64_overflow(long long base, long long exp, long long *res) {
122+
// special cases: n^0, 1^m, 0^m, -1^m, n^-m, NA^m, n^NA
123+
if (exp == 0) {
124+
*res = 1;
125+
return false;
126+
}
127+
if (base == 1) {
128+
*res = 1;
129+
return false;
130+
}
131+
if (base == 0) {
132+
if (exp < 0) {
133+
*res = NA_INTEGER64;
134+
return true;
135+
}
136+
*res = 0;
137+
return false;
138+
}
139+
if (base == -1) {
140+
*res = exp & 1 ? -1 : 1;
141+
return false;
142+
}
143+
if (exp < 0) {
144+
*res = 0;
145+
return false;
146+
}
147+
if (base == NA_INTEGER64 || exp == NA_INTEGER64) {
148+
*res = NA_INTEGER64;
149+
return false;
150+
}
151+
long long r = 1;
152+
while (exp > 0) {
153+
if (exp & 1 && mul64_overflow(r, base, &r)) {
154+
*res = NA_INTEGER64;
155+
return true;
156+
}
157+
exp >>= 1;
158+
if (exp > 0 && mul64_overflow(base, base, &base)) {
159+
*res = NA_INTEGER64;
160+
return true;
161+
}
162+
}
163+
*res = r;
164+
return false;
165+
}
166+
121167
#define PLUS64(e1,e2,ret,naflag) \
122168
if (e1 == NA_INTEGER64 || e2 == NA_INTEGER64) \
123169
ret = NA_INTEGER64; \
@@ -154,26 +200,6 @@ static inline bool mul64_overflow(long long a, long long b, long long *res) {
154200
ret = llroundl(longret); \
155201
}
156202

157-
#define POW64(e1,e2,ret,naflag,base,exp) \
158-
if (e1 == NA_INTEGER64 || e2 == NA_INTEGER64) \
159-
ret = NA_INTEGER64; \
160-
else if (e2 >= 0 || e1 == 1 || e1 == -1) { \
161-
ret = 1; \
162-
base = e1; \
163-
exp = e2; \
164-
while (exp > 0) { \
165-
if (exp % 2 == 1) { \
166-
if (mul64_overflow(ret, base, &ret)) { \
167-
ret = NA_INTEGER64; \
168-
naflag = TRUE; \
169-
break; \
170-
} \
171-
} \
172-
base *= base; \
173-
exp >>= 1; \
174-
} \
175-
}
176-
177203
#define POW64REAL(e1,e2,ret,naflag,longret) \
178204
if (e1 == NA_INTEGER64 || ISNAN(e2)) \
179205
ret = NA_INTEGER64; \

0 commit comments

Comments
 (0)