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