A068994

gistfile1.txt

#include "ulong_extras.h"
#include "nmod.h"
#include "fmpz.h"
#include "thread_support.h"

#define N_FAST1 UWORD(1000000)
#define N_FAST2 UWORD(100000000)
#define N_SLOW UWORD(10000000000000000000)

ulong * tab1;
ulong * tab2;

ulong
two_n_has_maybe_odd_digit(ulong e)
{
    fmpz_t m, x;
    ulong d;
    int res = 0;

    fmpz_init(x);
    fmpz_init(m);

    fmpz_ui_pow_ui(m, 10, 60);
    fmpz_set_ui(x, 2);
    fmpz_powm_ui(x, x, e, m);

    slong position = 0;

    while (!fmpz_is_zero(x))
    {
        d = fmpz_fdiv_ui(x, 10);

        if (d % 2)
        {
            if (position >= 37)
                flint_printf("%wu: position %wd\n", e, position);

            res = 1;
            break;
        }

        position++;

        fmpz_fdiv_q_ui(x, x, 10);
    }

    fmpz_clear(x);
    fmpz_clear(m);

    return res;
}

ulong
n_has_odd_digit(ulong x)
{
    ulong d;

    while (x != 0)
    {
        d = x % 10;

        if (d % 2)
            return 1;

        x = x / 10;
    }

    return 0;
}

void
compute_tab(ulong * tab, ulong N)
{
    ulong n;

    for (n = 0; n < N; n++)
        tab[n / FLINT_BITS] |= (n_has_odd_digit(n) << (n % FLINT_BITS));
}

ulong
check_range(ulong a, ulong b)
{
    ulong n, x1, x2, x3, count = 0;

    nmod_t n_fast1_mod;
    nmod_t n_fast2_mod;
    nmod_t n_slow_mod;

    nmod_init(&n_fast1_mod, N_FAST1);
    nmod_init(&n_fast2_mod, N_FAST2);
    nmod_init(&n_slow_mod, N_SLOW);

    x1 = nmod_pow_ui(2, a, n_fast1_mod);
    x2 = nmod_pow_ui(2, a, n_fast2_mod);

    for (n = a; n <= b; n++)
    {
        if (!(tab1[x1 / FLINT_BITS] & (UWORD(1) << (x1 % FLINT_BITS))))
        {
            if (!(tab2[x2 / FLINT_BITS] & (UWORD(1) << (x2 % FLINT_BITS))))
            {
                x3 = nmod_pow_ui(2, n, n_slow_mod);

                if (!n_has_odd_digit(x3))
                {
                    if (!two_n_has_maybe_odd_digit(n))
                    {
                        flint_printf("%wu: %wu\n", n, x3);
                        if (n > 11)
                            flint_abort();
                        count++;
                    }
                }
            }
        }

        x1 = n_addmod(x1, x1, N_FAST1);
        x2 = n_addmod(x2, x2, N_FAST2);
    }

    return count;
}

#define WORK_CHUNK_SIZE ((ulong) (1e9))

slong total_chunks = 0;
slong total_done = 0;

void
worker(slong i, void * args)
{
    ulong a, b;

    a = i * WORK_CHUNK_SIZE;
    b = (i + 1) * WORK_CHUNK_SIZE;

    ulong count = check_range(a, b);

    total_done++;
    flint_printf("      %wu in [%wu, %wu]    (~%.2f)\n", count, a, b, total_done / (double) total_chunks);
}

#define CHECK_UP_TO 1e13
#define NUM_THREADS 8

int main(void)
{
    tab1 = flint_calloc(sizeof(ulong), (N_FAST1 / FLINT_BITS + 1));
    tab2 = flint_calloc(sizeof(ulong), (N_FAST2 / FLINT_BITS + 1));

    compute_tab(tab1, N_FAST1);
    compute_tab(tab2, N_FAST2);

    total_chunks = (ulong) CHECK_UP_TO / WORK_CHUNK_SIZE;

    flint_set_num_threads(NUM_THREADS);
    flint_parallel_do(worker, NULL, (ulong) CHECK_UP_TO / WORK_CHUNK_SIZE, 0, 0);

    return 0;
}