This problem was posted on the GAP (Groups Algorithms & Programming) Forum some time ago. Roughly a week later, this partial solution was posted.

Knowing that GAP runs on an interpreted language on top of a C kernel, I thought I may be able to do better with C. After prototyping the situation in Python, my fears were realized: The digits associated with the orbit are HUGE. So, I decided to code in C using the GMP library and OpenMP. The main idea is that the orbit can be computed in parallel, going forwards and backwards until the orbit points agree. Not only will this cut the computation time roughly in half (when compared to a non-parallel C/GMP solution), but it should go much faster than GAP. And indeed, it does. My post to the forum announcing the first known solution can be found here. The maximum digit length found in the orbit is 76,785.

C/GMP/OpenMP Code

Here is the code. It requires the GNU MP Bignum Library (not in GCC) and OpenMP (in GCC). When the orbit points are within a specific digit length difference, only a single core continues the computation. Otherwise, both cores continue to compute the orbit in opposite directions. (There’s no makefile. I’ll give compiler instructions below.)

// Jason B. Hill
// Jason.B.Hill@Colorado.edu
// www.jasonbhill.com
 
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <omp.h>
 
//#define p 32
//#define p 736
//#define p 25952
#define p 173176
 
/*****************************************************************************/
/* Transpositions a,b,c and compositions g,g^-1                              */
/*****************************************************************************/
 
void a(mpz_t omega) {
    if(mpz_even_p(omega)==1) mpz_add_ui(omega, omega, 1);
    else mpz_sub_ui(omega, omega, 1);
}
 
void b(mpz_t omega) {
    if(mpz_congruent_ui_p(omega, 0, 5)==1) mpz_add_ui(omega, omega, 4);
    else if(mpz_congruent_ui_p(omega, 4, 5)==1) mpz_sub_ui(omega, omega, 4);
}
 
void c(mpz_t omega) {
    if(mpz_congruent_ui_p(omega, 1, 4)==1) {
        mpz_sub_ui(omega, omega, 1);
        mpz_divexact_ui(omega, omega, 4);
        mpz_mul_ui(omega, omega, 6);
    } else if(mpz_congruent_ui_p(omega, 0, 6)==1) {
        mpz_divexact_ui(omega, omega, 6);
        mpz_mul_ui(omega, omega, 4);
        mpz_add_ui(omega, omega, 1);
    }
}
 
void g(mpz_t omega) {
    a(omega);
    b(omega);
    c(omega);
}
 
void ginv(mpz_t omega) {
    c(omega);
    b(omega);
    a(omega);
}
 
/*****************************************************************************/
/* Main                                                                      */
/*****************************************************************************/
 
int main(void) {
    unsigned long       n = p;
    unsigned long long  i0 = 0;
    unsigned long long  i1 = 0;
    int                 th_id;
    _Bool               sstop = 0;
    size_t              s = 0;
    size_t              c0, c1;
    mpz_t               omega0, omega1;
 
    omp_set_num_threads(2);
 
    mpz_init(omega0);
    mpz_init(omega1);
 
    mpz_set_ui(omega0, n);
    mpz_set_ui(omega1, n);
 
    c0 = mpz_sizeinbase(omega0, 10);
    c1 = mpz_sizeinbase(omega1, 10);
 
    #pragma omp parallel private(th_id) \
    shared(omega0,omega1,s,sstop,c0,c1,i0,i1)
    {
        th_id = omp_get_thread_num();
 
        if(th_id == 1) {
            g(omega1);
            i1++;
            if(mpz_cmp(omega0,omega1)==0) sstop = 1;
        }
 
        #pragma omp barrier
 
        while(!sstop) {
            if(th_id == 0) {
                if(abs(c0 - c1) > 20) {
                    ginv(omega0);
                    c0 = mpz_sizeinbase(omega0, 10);
                    i0++;
                }
            } else if(th_id == 1) {
                if(abs(c0 - c1) > 5) {
                    g(omega1);
                    c1 = mpz_sizeinbase(omega1, 10);
                    i1++;
                } else {
                    if(mpz_cmp(omega0,omega1)==0) sstop = 1;
                    else {
                        g(omega1);
                        c1 = mpz_sizeinbase(omega1, 10);
                        i1++;
                    }
                }
            }
 
            #pragma omp flush(sstop,c0,c1,i0,i1)
 
            if(th_id == 0) {
                if(c0 > s) {
                    s = c0;
                    printf("Core 0: digit length increased to %ld\n", s);
                    printf("iterations: %lld (core0) %lld (core1) %lld (total)\
\n\n",i0,i1,i0+i1);
                }
                if(c1 > s) {
                    s = c1;
                    printf("Core 1: digit length increased to %ld\n", s);
                    printf("iterations: %lld (core0) %lld (core1) %lld (total)\
\n\n",i0,i1,i0+i1);
                }
                if((i0+i1)%100000000==0) {
                    printf("digit length: %ld\n", s);
                    printf("iterations: %lld (core0) %lld (core1) %lld (total)\
\n\n",i0,i1,i0+i1);
                }
            }
        }
    }
 
    printf("total iterations: %lld\n",i0+i1);
    mpz_clear(omega0);
    mpz_clear(omega1);
 
    return 0;
}

Using GCC with OpenMP and GMP, one can compile the code on a multicore machine with, for example, the following command.

gcc -O3 -o length-residue-orbit-omp length-residue-orbit-omp.c -lgmp -fopenmp

Success!

The result: total iterations: 47610700792, is returned in roughly 3.1 days on a 2.9 GHz 3rd generation Intel Core i7 (i7-3520M).

Posted in C, OpenMP
Share this post, let the world know

Leave a Reply

Your email address will not be published. Required fields are marked *

*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> <pre lang="" line="" escaped="" highlight="">