/* Copyright (C) 2016 Pascal Molin 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 "dlog.h" #include #if FLINT_BITS == 64 #define LIM UWORD(1000000000000) #else #define LIM UWORD(1000000000) #endif int main() { slong iter; flint_rand_t state; flint_printf("modpe...."); fflush(stdout); flint_randinit(state); for (iter = 0; iter < 1000; iter++) { ulong p, e, pe, a; if (iter == 0) { p = 2; pe = 8; a = 5; e = 3; } else { p = pe = n_randprime(state, 10, 0); a = (p == 40487) ? 10 : n_primitive_root_prime(p); e = 1; } for (; pe < LIM; pe *= p, e++) { ulong k, phi; nmod_t mod; dlog_modpe_t modpe; nmod_init(&mod, pe); phi = (p == 2) ? pe / 4 : pe - pe / p; dlog_modpe_init(modpe, a, p, e, pe, 10); for (k = 0; k < 100 && k < p; k++) { ulong l, b, x; l = n_randint(state, phi); b = nmod_pow_ui(a, l, mod); if ((x = dlog_modpe(modpe, b)) != l) { flint_printf("FAIL modpe: %wu^%wu = %wu [%wu^%wu]\n\n", a, l, b, p, e); flint_printf("modpe returned %wu\n\n", x); flint_abort(); } } dlog_modpe_clear(modpe); /* multiplication can overflow on 32-bit */ if ((double) pe * p > LIM) break; } } flint_randclear(state); flint_cleanup(); flint_printf("PASS\n"); return EXIT_SUCCESS; }