Skip to content

Commit 497769e

Browse files
shibatchshibatch
shibatch
and
shibatch
authored
[Quad] Add frexp and modf
This patch adds quad-precision frexp and modf. Co-authored-by: shibatch <[email protected]>
1 parent f9d5b33 commit 497769e

File tree

8 files changed

+312
-13
lines changed

8 files changed

+312
-13
lines changed

‎src/common/keywords.txt

+1
Original file line numberDiff line numberDiff line change
@@ -674,3 +674,4 @@ rint_tdx_tdx
674674
fmod_tdx_tdx_tdx
675675
remainder_tdx_tdx_tdx
676676
cbrt_tdx_tdx
677+
frexp_tdx_tdx_pvi

‎src/quad-tester/qiutcuda.cu

+36
Original file line numberDiff line numberDiff line change
@@ -98,6 +98,9 @@ __global__ void xfdimq_u05(Sleef_quadx1 *r, Sleef_quadx1 *a0, Sleef_quadx1 *a1)
9898
__global__ void xfmodq(Sleef_quadx1 *r, Sleef_quadx1 *a0, Sleef_quadx1 *a1) { *r = Sleef_fmodq1_cuda(*a0, *a1); }
9999
__global__ void xremainderq(Sleef_quadx1 *r, Sleef_quadx1 *a0, Sleef_quadx1 *a1) { *r = Sleef_remainderq1_cuda(*a0, *a1); }
100100

101+
__global__ void xfrexpq(Sleef_quadx1 *r, Sleef_quadx1 *a0, int *i0) { *r = Sleef_frexpq1_cuda(*a0, i0); }
102+
__global__ void xmodfq(Sleef_quadx1 *r, Sleef_quadx1 *a0, Sleef_quadx1 *a1) { *r = Sleef_modfq1_cuda(*a0, a1); }
103+
101104
__global__ void xtruncq(Sleef_quadx1 *r, Sleef_quadx1 *a0) { *r = Sleef_truncq1_cuda(*a0); }
102105
__global__ void xfloorq(Sleef_quadx1 *r, Sleef_quadx1 *a0) { *r = Sleef_floorq1_cuda(*a0); }
103106
__global__ void xceilq(Sleef_quadx1 *r, Sleef_quadx1 *a0) { *r = Sleef_ceilq1_cuda(*a0); }
@@ -247,6 +250,37 @@ typedef union {
247250
} \
248251
}
249252

253+
#define func_q_q_pi(funcStr, funcName) { \
254+
while (startsWith(buf, funcStr " ")) { \
255+
sentinel = 0; \
256+
cnv128 c0; \
257+
sscanf(buf, funcStr " %" PRIx64 ":%" PRIx64, &c0.h, &c0.l); \
258+
*a0 = Sleef_setq1_cuda(*a0, 0, c0.q); \
259+
funcName<<<1, 1>>>(r, a0, i0); \
260+
cudaDeviceSynchronize(); \
261+
c0.q = Sleef_getq1_cuda(*r, 0); \
262+
printf("%" PRIx64 ":%" PRIx64 " %d\n", c0.h, c0.l, *i0); \
263+
fflush(stdout); \
264+
if (fgets(buf, BUFSIZE-1, stdin) == NULL) break; \
265+
} \
266+
}
267+
268+
#define func_q_q_pq(funcStr, funcName) { \
269+
while (startsWith(buf, funcStr " ")) { \
270+
sentinel = 0; \
271+
cnv128 c0, c1; \
272+
sscanf(buf, funcStr " %" PRIx64 ":%" PRIx64, &c0.h, &c0.l); \
273+
*a0 = Sleef_setq1_cuda(*a0, 0, c0.q); \
274+
funcName<<<1, 1>>>(r, a0, a1); \
275+
cudaDeviceSynchronize(); \
276+
c0.q = Sleef_getq1_cuda(*r, 0); \
277+
c1.q = Sleef_getq1_cuda(*a1, 0); \
278+
printf("%" PRIx64 ":%" PRIx64 " %" PRIx64 ":%" PRIx64 "\n", c0.h, c0.l, c1.h, c1.l); \
279+
fflush(stdout); \
280+
if (fgets(buf, BUFSIZE-1, stdin) == NULL) break; \
281+
} \
282+
}
283+
250284
int main(int argc, char **argv) {
251285
#if 0
252286
cuInit(0);
@@ -327,6 +361,8 @@ int main(int argc, char **argv) {
327361
func_q_q_q("fdimq_u05", xfdimq_u05);
328362
func_q_q_q("fmodq", xfmodq);
329363
func_q_q_q("remainderq", xremainderq);
364+
func_q_q_pi("frexpq", xfrexpq);
365+
func_q_q_pq("modfq", xmodfq);
330366

331367
func_q_q("truncq", xtruncq);
332368
func_q_q("floorq", xfloorq);

‎src/quad-tester/qiutsimd.c

+40
Original file line numberDiff line numberDiff line change
@@ -406,6 +406,44 @@ typedef union {
406406
} \
407407
}
408408

409+
#define func_q_q_pi(funcStr, funcName) { \
410+
while (startsWith(buf, funcStr " ")) { \
411+
sentinel = 0; \
412+
int lane = xrand() % VECTLENDP; \
413+
cnv128 c0; \
414+
sscanf(buf, funcStr " %" PRIx64 ":%" PRIx64, &c0.h, &c0.l); \
415+
VARGQUAD a0; \
416+
memrand(&a0, SIZEOF_VARGQUAD); \
417+
a0 = xsetq(a0, lane, c0.q); \
418+
vint vi; \
419+
a0 = funcName(a0, &vi); \
420+
c0.q = xgetq(a0, lane); \
421+
int t[VECTLENDP*2]; \
422+
vstoreu_v_p_vi(t, vi); \
423+
printf("%" PRIx64 ":%" PRIx64 " %d\n", c0.h, c0.l, t[lane]); \
424+
fflush(stdout); \
425+
if (fgets(buf, BUFSIZE-1, stdin) == NULL) break; \
426+
} \
427+
}
428+
429+
#define func_q_q_pq(funcStr, funcName) { \
430+
while (startsWith(buf, funcStr " ")) { \
431+
sentinel = 0; \
432+
int lane = xrand() % VECTLENDP; \
433+
cnv128 c0, c1; \
434+
sscanf(buf, funcStr " %" PRIx64 ":%" PRIx64, &c0.h, &c0.l); \
435+
VARGQUAD a0, a1; \
436+
memrand(&a0, SIZEOF_VARGQUAD); \
437+
a0 = xsetq(a0, lane, c0.q); \
438+
a0 = funcName(a0, &a1); \
439+
c0.q = xgetq(a0, lane); \
440+
c1.q = xgetq(a1, lane); \
441+
printf("%" PRIx64 ":%" PRIx64 " %" PRIx64 ":%" PRIx64 "\n", c0.h, c0.l, c1.h, c1.l); \
442+
fflush(stdout); \
443+
if (fgets(buf, BUFSIZE-1, stdin) == NULL) break; \
444+
} \
445+
}
446+
409447
#define func_strtoq(funcStr) { \
410448
while (startsWith(buf, funcStr " ")) { \
411449
sentinel = 0; \
@@ -562,6 +600,8 @@ int do_test(int argc, char **argv) {
562600
func_q_q_q("fdimq_u05", xfdimq_u05);
563601
func_q_q_q("fmodq", xfmodq);
564602
func_q_q_q("remainderq", xremainderq);
603+
func_q_q_pi("frexpq", xfrexpq);
604+
func_q_q_pq("modfq", xmodfq);
565605

566606
func_q_q("truncq", xtruncq);
567607
func_q_q("floorq", xfloorq);

‎src/quad-tester/qtester.c

+164
Original file line numberDiff line numberDiff line change
@@ -198,6 +198,31 @@ typedef union {
198198
return c.q; \
199199
} while(0)
200200

201+
#define child_q_q_pi(funcStr, arg) do { \
202+
char str[256]; \
203+
cnv128 c; \
204+
c.q = arg; \
205+
sprintf(str, funcStr " %" PRIx64 ":%" PRIx64 "\n", c.h, c.l); \
206+
write(ptoc[1], str, strlen(str)); \
207+
if (fgets(str, 255, fpctop) == NULL) stop("child " funcStr); \
208+
int i; \
209+
sscanf(str, "%" PRIx64 ":%" PRIx64 " %d", &c.h, &c.l, &i); \
210+
*ptr = i; \
211+
return c.q; \
212+
} while(0)
213+
214+
#define child_q_q_pq(funcStr, arg) do { \
215+
char str[256]; \
216+
cnv128 c0, c1; \
217+
c0.q = arg; \
218+
sprintf(str, funcStr " %" PRIx64 ":%" PRIx64 "\n", c0.h, c0.l); \
219+
write(ptoc[1], str, strlen(str)); \
220+
if (fgets(str, 255, fpctop) == NULL) stop("child " funcStr); \
221+
sscanf(str, "%" PRIx64 ":%" PRIx64 " %" PRIx64 ":%" PRIx64, &c0.h, &c0.l, &c1.h, &c1.l); \
222+
*ptr = c1.q; \
223+
return c0.q; \
224+
} while(0)
225+
201226
#define child_q_str(funcStr, arg) do { \
202227
char str[256]; \
203228
sprintf(str, funcStr " %s\n", arg); \
@@ -283,6 +308,9 @@ Sleef_quad child_ceilq(Sleef_quad x) { child_q_q("ceilq", x); }
283308
Sleef_quad child_roundq(Sleef_quad x) { child_q_q("roundq", x); }
284309
Sleef_quad child_rintq(Sleef_quad x) { child_q_q("rintq", x); }
285310

311+
Sleef_quad child_frexpq(Sleef_quad x, int *ptr) { child_q_q_pi("frexpq", x); }
312+
Sleef_quad child_modfq(Sleef_quad x, Sleef_quad *ptr) { child_q_q_pq("modfq", x); }
313+
286314
//
287315

288316
#define cmpDenorm_q(mpfrFunc, childFunc, argx) do { \
@@ -326,6 +354,36 @@ Sleef_quad child_rintq(Sleef_quad x) { child_q_q("rintq", x); }
326354
} \
327355
} while(0)
328356

357+
#define cmpDenorm_q_pi(mpfrFunc, childFunc, argx) do { \
358+
mpfr_set_f128(frx, argx, GMP_RNDN); \
359+
mpfr_exp_t e; \
360+
mpfrFunc(&e, frz, frx, GMP_RNDN); \
361+
int i; \
362+
Sleef_quad t = childFunc(argx, &i); \
363+
double u = countULPf128(t, frz, 1); \
364+
if (u >= 10 || i != (int)e) { \
365+
fprintf(stderr, "\narg = %s\ntest = %s, %d\ncorrect = %s, %d\nulp = %g\n", \
366+
sprintf128(argx), sprintf128(t), i, sprintfr(frz), (int)e, u); \
367+
success = 0; \
368+
break; \
369+
} \
370+
} while(0)
371+
372+
#define cmpDenorm_q_pq(mpfrFunc, childFunc, argx) do { \
373+
mpfr_set_f128(frx, argx, GMP_RNDN); \
374+
mpfrFunc(fry, frz, frx, GMP_RNDN); \
375+
Sleef_quad qi, qf; \
376+
qf = childFunc(argx, &qi); \
377+
double u = countULPf128(qf, frz, 1); \
378+
double v = countULPf128(qi, fry, 1); \
379+
if (u >= 10 || v >= 10) { \
380+
fprintf(stderr, "\narg = %s\ntest = %s, %s\ncorrect = %s, %s\nulp = %g, %g\n", \
381+
sprintf128(argx), sprintf128(qf), sprintf128(qi), sprintfr(frz), sprintfr(fry), u, v); \
382+
success = 0; \
383+
break; \
384+
} \
385+
} while(0)
386+
329387
#define checkAccuracy_q(mpfrFunc, childFunc, argx, bound) do { \
330388
mpfr_set_f128(frx, argx, GMP_RNDN); \
331389
mpfrFunc(frz, frx, GMP_RNDN); \
@@ -369,6 +427,39 @@ Sleef_quad child_rintq(Sleef_quad x) { child_q_q("rintq", x); }
369427
} \
370428
} while(0)
371429

430+
#define checkAccuracy_q_pi(mpfrFunc, childFunc, argx, bound) do { \
431+
mpfr_set_f128(frx, argx, GMP_RNDN); \
432+
mpfr_exp_t ex; \
433+
mpfrFunc(&ex, frz, frx, GMP_RNDN); \
434+
int i; \
435+
Sleef_quad t = childFunc(argx, &i); \
436+
double e = countULPf128(t, frz, 0); \
437+
maxError = fmax(maxError, e); \
438+
if (e > bound || i != (int)ex) { \
439+
fprintf(stderr, "\narg = %s, test = %s, %d, correct = %s, %d, ULP = %lf\n", \
440+
sprintf128(argx), sprintf128(t), i, sprintfr(frz), (int)ex, countULPf128(t, frz, 0)); \
441+
success = 0; \
442+
break; \
443+
} \
444+
} while(0)
445+
446+
#define checkAccuracy_q_pq(mpfrFunc, childFunc, argx, bound) do { \
447+
mpfr_set_f128(frx, argx, GMP_RNDN); \
448+
mpfrFunc(fry, frz, frx, GMP_RNDN); \
449+
Sleef_quad qi, qf; \
450+
qf = childFunc(argx, &qi); \
451+
double ef = countULPf128(qf, frz, 0); \
452+
double ei = countULPf128(qi, fry, 0); \
453+
maxError = fmax(maxError, ef); \
454+
maxError = fmax(maxError, ei); \
455+
if (ef > bound || ei > bound) { \
456+
fprintf(stderr, "\narg = %s, test = %s, %s, correct = %s, %s, ULP = %lf, %lf\n", \
457+
sprintf128(argx), sprintf128(qf), sprintf128(qi), sprintfr(frz), sprintfr(fry), ef, ei); \
458+
success = 0; \
459+
break; \
460+
} \
461+
} while(0)
462+
372463
#define testComparison(mpfrFunc, childFunc, argx, argy) do { \
373464
mpfr_set_f128(frx, argx, GMP_RNDN); \
374465
mpfr_set_f128(fry, argy, GMP_RNDN); \
@@ -408,6 +499,20 @@ Sleef_quad child_rintq(Sleef_quad x) { child_q_q("rintq", x); }
408499
} \
409500
} while(0)
410501

502+
#define cmpDenormOuterLoop_q_pi(mpfrFunc, childFunc, checkVals) do { \
503+
for(int i=0;i<sizeof(checkVals)/sizeof(char *);i++) { \
504+
Sleef_quad a0 = cast_q_str(checkVals[i]); \
505+
cmpDenorm_q_pi(mpfrFunc, childFunc, a0); \
506+
} \
507+
} while(0)
508+
509+
#define cmpDenormOuterLoop_q_pq(mpfrFunc, childFunc, checkVals) do { \
510+
for(int i=0;i<sizeof(checkVals)/sizeof(char *);i++) { \
511+
Sleef_quad a0 = cast_q_str(checkVals[i]); \
512+
cmpDenorm_q_pq(mpfrFunc, childFunc, a0); \
513+
} \
514+
} while(0)
515+
411516
//
412517

413518
#define checkAccuracyOuterLoop_q_q(mpfrFunc, childFunc, minStr, maxStr, sign, nLoop, bound, seed) do { \
@@ -471,6 +576,38 @@ Sleef_quad child_rintq(Sleef_quad x) { child_q_q("rintq", x); }
471576
} \
472577
} while(0)
473578

579+
#define checkAccuracyOuterLoop_q_pi(mpfrFunc, childFunc, minStr, maxStr, sign, nLoop, bound, seed) do { \
580+
xsrand(seed); \
581+
Sleef_quad min = cast_q_str(minStr), max = cast_q_str(maxStr); \
582+
for(int i=0;i<nLoop && success;i++) { \
583+
Sleef_quad x = rndf128(min, max, sign); \
584+
checkAccuracy_q_pi(mpfrFunc, childFunc, x, bound); \
585+
} \
586+
} while(0)
587+
588+
#define checkAccuracyOuterLoop2_q_pi(mpfrFunc, childFunc, checkVals, bound) do { \
589+
for(int i=0;i<sizeof(checkVals)/sizeof(char *);i++) { \
590+
Sleef_quad x = cast_q_str(checkVals[i]); \
591+
checkAccuracy_q_pi(mpfrFunc, childFunc, x, bound); \
592+
} \
593+
} while(0)
594+
595+
#define checkAccuracyOuterLoop_q_pq(mpfrFunc, childFunc, minStr, maxStr, sign, nLoop, bound, seed) do { \
596+
xsrand(seed); \
597+
Sleef_quad min = cast_q_str(minStr), max = cast_q_str(maxStr); \
598+
for(int i=0;i<nLoop && success;i++) { \
599+
Sleef_quad x = rndf128(min, max, sign); \
600+
checkAccuracy_q_pq(mpfrFunc, childFunc, x, bound); \
601+
} \
602+
} while(0)
603+
604+
#define checkAccuracyOuterLoop2_q_pq(mpfrFunc, childFunc, checkVals, bound) do { \
605+
for(int i=0;i<sizeof(checkVals)/sizeof(char *);i++) { \
606+
Sleef_quad x = cast_q_str(checkVals[i]); \
607+
checkAccuracy_q_pq(mpfrFunc, childFunc, x, bound); \
608+
} \
609+
} while(0)
610+
474611
//
475612

476613
void checkResult(int success, double e) {
@@ -545,6 +682,17 @@ void do_test(int options) {
545682
"NaN"
546683
};
547684

685+
static const char *finiteCheckVals[] = {
686+
"-0.0", "0.0", "+0.25", "-0.25", "+0.5", "-0.5", "+0.75", "-0.75", "+1.0", "-1.0",
687+
"+1.25", "-1.25", "+1.5", "-1.5", "+2.0", "-2.0", "+2.5", "-2.5", "+3.0", "-3.0",
688+
"+4.0", "-4.0", "+5.0", "-5.0", "+6.0", "-6.0", "+7.0", "-7.0",
689+
"1.234", "-1.234", "+1.234e+100", "-1.234e+100", "+1.234e-100", "-1.234e-100",
690+
"+1.234e+3000", "-1.234e+3000", "+1.234e-3000", "-1.234e-3000",
691+
"3.1415926535897932384626433832795028841971693993751058209749445923078164",
692+
"+" STR_QUAD_MIN, "-" STR_QUAD_MIN,
693+
"+" STR_QUAD_DENORM_MIN, "-" STR_QUAD_DENORM_MIN,
694+
};
695+
548696
static const char *trigCheckVals[] = {
549697
"3.141592653589793238462643383279502884197169399375105820974944592307",
550698
"6.283185307179586476925286766559005768394338798750211641949889184615",
@@ -665,6 +813,22 @@ void do_test(int options) {
665813

666814
//
667815

816+
fprintf(stderr, "frexp : ");
817+
maxError = 0;
818+
cmpDenormOuterLoop_q_pi(mpfr_frexp, child_frexpq, finiteCheckVals);
819+
checkAccuracyOuterLoop2_q_pi(mpfr_frexp, child_frexpq, finiteCheckVals, 0);
820+
checkAccuracyOuterLoop_q_pi(mpfr_frexp, child_frexpq, "1e-4000", "1e+4000", 1, 10 * NTEST, 0, 1);
821+
checkResult(success, maxError);
822+
823+
fprintf(stderr, "modf : ");
824+
maxError = 0;
825+
cmpDenormOuterLoop_q_pq(mpfr_modf, child_modfq, stdCheckVals);
826+
checkAccuracyOuterLoop2_q_pq(mpfr_modf, child_modfq, stdCheckVals, 0);
827+
checkAccuracyOuterLoop_q_pq(mpfr_modf, child_modfq, "1e-4000", "1e+4000", 1, 10 * NTEST, 0, 1);
828+
checkResult(success, maxError);
829+
830+
//
831+
668832
fprintf(stderr, "cast_from_doubleq : ");
669833
{
670834
xsrand(0);

‎src/quad-tester/tester2simdqp.c

+16-1
Original file line numberDiff line numberDiff line change
@@ -441,7 +441,22 @@ int main(int argc,char **argv)
441441
}
442442
}
443443

444-
if (cnt % 7 == 0) {
444+
{
445+
mpfr_modf(frz, frw, frx, GMP_RNDN);
446+
447+
a2 = xmodfq(a0, &a3);
448+
double u0 = countULPf128(q2 = vget(a2, e), frw, 0);
449+
double u1 = countULPf128(q3 = vget(a3, e), frz, 0);
450+
451+
if (u0 > 0 || u1 > 0) {
452+
printf(ISANAME " modf arg=%s ulp=%.20g, %.20g\n", sprintf128(q0), u0, u1);
453+
printf("test = %s, %s\n", sprintf128(q2), sprintf128(q3));
454+
printf("corr = %s, %s\n\n", sprintf128(mpfr_get_f128(frw, GMP_RNDN)), sprintf128(mpfr_get_f128(frz, GMP_RNDN)));
455+
fflush(stdout); ecnt++;
456+
}
457+
}
458+
459+
if (cnt % 101 == 0) {
445460
{
446461
mpfr_fmod(frz, frx, fry, GMP_RNDN);
447462

0 commit comments

Comments
 (0)