/* Copyright (C) 2012 Fredrik Johansson This file is part of Arb. Arb is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. See . */ #include "arf.h" int arf_mul_naive(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd) { if (rnd == ARF_RND_NEAR) { arf_mul(z, x, y, ARF_PREC_EXACT, rnd); return arf_set_round(z, z, prec, rnd); } else { fmpr_t a, b; slong r; fmpr_init(a); fmpr_init(b); arf_get_fmpr(a, x); arf_get_fmpr(b, y); r = fmpr_mul_naive(a, a, b, prec, rnd); arf_set_fmpr(z, a); fmpr_clear(a); fmpr_clear(b); return (r == FMPR_RESULT_EXACT) ? 0 : 1; } } int main() { slong iter, iter2; flint_rand_t state; flint_printf("mul_via_mpfr...."); fflush(stdout); flint_randinit(state); for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++) { arf_t x, y, z, v; slong prec, r1, r2; arf_rnd_t rnd; arf_init(x); arf_init(y); arf_init(z); arf_init(v); for (iter2 = 0; iter2 < 100; iter2++) { arf_randtest_special(x, state, 2000, 100); arf_randtest_special(y, state, 2000, 100); prec = 2 + n_randint(state, 2000); if (n_randint(state, 50) == 0) prec = ARF_PREC_EXACT; switch (n_randint(state, 5)) { case 0: rnd = ARF_RND_DOWN; break; case 1: rnd = ARF_RND_UP; break; case 2: rnd = ARF_RND_FLOOR; break; case 3: rnd = ARF_RND_CEIL; break; default: rnd = ARF_RND_NEAR; break; } switch (n_randint(state, 5)) { case 0: r1 = arf_mul_via_mpfr(z, x, y, prec, rnd); r2 = arf_mul_naive(v, x, y, prec, rnd); if (!arf_equal(z, v) || r1 != r2) { flint_printf("FAIL!\n"); flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); flint_printf("x = "); arf_print(x); flint_printf("\n\n"); flint_printf("y = "); arf_print(y); flint_printf("\n\n"); flint_printf("z = "); arf_print(z); flint_printf("\n\n"); flint_printf("v = "); arf_print(v); flint_printf("\n\n"); flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); flint_abort(); } break; case 1: r1 = arf_mul_via_mpfr(z, x, x, prec, rnd); r2 = arf_mul_naive(v, x, x, prec, rnd); if (!arf_equal(z, v) || r1 != r2) { flint_printf("FAIL (aliasing 1)!\n"); flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); flint_printf("x = "); arf_print(x); flint_printf("\n\n"); flint_printf("z = "); arf_print(z); flint_printf("\n\n"); flint_printf("v = "); arf_print(v); flint_printf("\n\n"); flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); flint_abort(); } break; case 2: r2 = arf_mul_naive(v, x, x, prec, rnd); r1 = arf_mul_via_mpfr(x, x, x, prec, rnd); if (!arf_equal(v, x) || r1 != r2) { flint_printf("FAIL (aliasing 2)!\n"); flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); flint_printf("x = "); arf_print(x); flint_printf("\n\n"); flint_printf("z = "); arf_print(z); flint_printf("\n\n"); flint_printf("v = "); arf_print(v); flint_printf("\n\n"); flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); flint_abort(); } break; case 3: r2 = arf_mul_naive(v, x, y, prec, rnd); r1 = arf_mul_via_mpfr(x, x, y, prec, rnd); if (!arf_equal(x, v) || r1 != r2) { flint_printf("FAIL (aliasing 3)!\n"); flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); flint_printf("x = "); arf_print(x); flint_printf("\n\n"); flint_printf("y = "); arf_print(y); flint_printf("\n\n"); flint_printf("v = "); arf_print(v); flint_printf("\n\n"); flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); flint_abort(); } break; default: r2 = arf_mul_naive(v, x, y, prec, rnd); r1 = arf_mul_via_mpfr(x, y, x, prec, rnd); if (!arf_equal(x, v) || r1 != r2) { flint_printf("FAIL (aliasing 4)!\n"); flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); flint_printf("x = "); arf_print(x); flint_printf("\n\n"); flint_printf("y = "); arf_print(y); flint_printf("\n\n"); flint_printf("v = "); arf_print(v); flint_printf("\n\n"); flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); flint_abort(); } break; } } arf_clear(x); arf_clear(y); arf_clear(z); arf_clear(v); } flint_randclear(state); flint_cleanup(); flint_printf("PASS\n"); return EXIT_SUCCESS; }