A000057, A064414 and more.

[Daniel Mondot]

Hi everyone.

I recently became interested in A000057 (primes dividing all Fibonacci sequences), and tried to improve the b-file. My first attempt was really slow, as I was doing a test that was taking a time proportional to k^2. Then I realized that the program by Charles Greathouse, which I rewrote in C because I don't speak PARI, was really fast since the test takes a time that is proportional to k instead of k^2. So it was easy to boost the b-file to 10000 terms.

Then I realized that A000057 is just a subsequence of A064414, which should probably be renamed just "Numbers dividing all Fibonacci sequences". However the test that worked for A000057 doesn't work for A064414. The composite numbers need a test that takes a time proportional to k^2. However, to get more terms, we can be smart about it:

- We can separate primes from composite numbers and use the fast test on the primes, and the slow test on the composite numbers.

- If a number doesn't divide all fibonacci sequences, then neither do any of its multiples. So, as we are done processing each number, we can update a sieve to eliminate all the future numbers that we know can't be in the sequence.

- and lastly, except for numbers that are being eliminated in a sieve, each new number in the sequence doesn't depend on previous numbers. So we take advantage of today's newer computers having multiple cores, and have each core process a different number in parallel.

So I have been experimenting with multi-threaded programming. After many trial and errors, I was able to get 10000 terms for A064414 on a 4-core computer in 7.5 days. Then I reran a slightly modified version of this program on a 32-core computer and got the same result in less than 8 hours (about 24 times faster). 

Of course my program is much more complicated than a single-thread version. Not only the next number needs to be assigned (without duplication) to each thread as they finish their previous task, but the time it takes for each test to complete varies wildly, so a single thread is responsible to collect all the results, and is responsible to display them in the right order. I may be wrong, but I don't believe that programs written in PARI or Mathematica could take advantage of multithreading.

Now that I have 10000 terms for both A000057 and A064414, I can do some analysis.

We already know that except for 2, all terms of A000057 are congruent to 3 or 7 modulo 20, but what about A064414?

It turns out that there are some interesting patterns, for instance:

- if you look at the terms modulo 20, most terms are congruent to 3,6,7 and 14. Exceptions are at 1,2,4, 9, and 18. None of them are congruent to 5, 10,11,12,13,15,16,17 or 19.

- if you look at the terms modulo 24, most terms are congruent to 7,11,14, 19 and 23. Exceptions are at 1,2,3, 4, 6, and 9. None of them are congruent to 5, 8, 10, 12, 13, 15, 16, 17, 18, 20, 21 or 22.

Also, almost none of them are multiple of 3 or7.

Now if you look at the composite numbers of A064414 only, the sequence begins 1, 4, 6, 9, 14, 27, 49, 81, 86, 98, etc... and contains a lot of powers.

Then take the first difference of these numbers modulo 24, and you get this:

You can probably see a horizontal symmetry around 12, and a solid line at the bottom.

That means that:

- Most numbers are multiple of 24.

- The exceptions are relatively rare.

- None of them are congruent to 1, 4, 6, 10, 14, 16, 18, 20, 21, or 23 mod 24.

- And with the exceptions of 4, 9, 27 and 81, all exceptions are located between 2 multiples of 24, there appears to be no double exceptions back to back.

Now, I currently only have about 3000 terms, so there might still be more exceptions further down the sequence.

We can also look at the first difference of the composite terms of A064414 modulo other values. For instance, the graph for modulo 17 is curious:

Some values seem a lot more rare than others, for no obvious reasons.

Similar graphs for 41, 47, etc...

General questions: 
Would it be worth having 5 new sequences:

1) The composite numbers of A064414 (I currently have about 3000 terms, hoping to get 10000 by the end of the week).

2) The first differences of A064414,

3) The first differences of its prime terms (A000057),

4) The first difference of its composite terms (new sequence (1)) ?

5) And since about 3/4 of the terms of A064414 are prime, and most of what remains are multiple of 24, we could have a sequence with just the remaining exceptions: the terms of A064414 that are neither primes or multiple of 24.

PS: I have uploaded new 10000 terms b-files for A000057 and A064414, but have not submitted them in case I need to make some additional changes to the comments.

Daniel.

[Giovanni Resta]

Hi Daniel,

since you are interested in these sequences, you may find useful the one I wrote for A00057.
It is a very rudimental  short C program, but computes the terms < 10^9 in about 8 seconds on my 10+ years old 6 cores CPU.

It uses Kim Walisch primesieve library for prime generation (it is a bit overkill for such small numbers).

It uses OpenMP for parallelism.

Giovanni

#include<stdio.h>
#include<stdlib.h>
#include<inttypes.h>
#include<omp.h>
#include<primesieve.h>

// NOTE: ONLY FOR NUMBERS < 10^9
// primesieve library from http://primesieve.org/ (Kim Walisch)
// compiled under linux with
// gcc -O3 a57.c -fopenmp -lprimesieve -o a57

// for factorizing numbers up to 10^9
// note that Primes[NP-1]^2 > 10^9
#define NP 3500
uint32_t  Primes[NP];

// computes A.B  mod m
inline void Mul(uint32_t A[2][2], uint32_t B[2][2], uint32_t m){
  uint64_t t11,t12,t21,t22;
  t11 = A[0][0]*(uint64_t)B[0][0]+A[0][1]*(uint64_t)B[1][0];
  t12 = A[0][0]*(uint64_t )B[0][1]+A[0][1]*(uint64_t)B[1][1];
  t21 = A[1][0]*(uint64_t )B[0][0]+A[1][1]*(uint64_t)B[1][0];
  t22 = A[1][0]*(uint64_t )B[0][1]+A[1][1]*(uint64_t)B[1][1];
  A[0][0]=t11%m; A[0][1]=t12%m; A[1][0]=t21%m; A[1][1]=t22%m;
}

// computes Fibonacci(e) modulo m for uint32_t values

// using matrix exponentiation via repeated squaring
uint32_t fimod(uint32_t e, uint32_t m){
  uint32_t X[2][2]={{1,0},{0,1}};
  uint32_t Y[2][2]={{1,1},{1,0}};
  while (e) {
      if (e&1) Mul(X,Y,m);
      Mul(Y,Y,m);
      e>>=1;
    }
  return X[0][1];
}

uint32_t Check(uint32_t n){
  uint32_t m=n+1;

 // check factor 2
  if ((m&1)==0 && fimod(m/2,n)==0) return 0;
  // makes m odd
  while ((m&1)==0) m >>= 1;
  // trivial factorization, no need to store exponents
  int j=1;
  while (Primes[j]*Primes[j] <= m){    
    if (0==(m % Primes[j])){
      if (fimod((n+1)/Primes[j],n)==0)
return 0;
      // delete prime factors
      do m /= Primes[j]; while (0==(m%Primes[j]));
    }
    j++;
  }
  // last prime factor
  if (m > 1 && fimod((n+1)/m,n)==0) return 0;
  // else...
  return 1;
}  

// I processs this many primes in parallel
#define BLOCK (10000)

int main(){
 uint32_t Candidate[BLOCK]; // to be checked
 uint32_t Res[BLOCK]; // result of checking
 
 primesieve_iterator it;
 primesieve_init(&it);
 // stores primes for factorization
 for (int i=0; i < NP; i++)
   Primes[i] = primesieve_next_prime(&it); 
 // set back prime sieve
 primesieve_jump_to(&it, 0, 1.2e9);
 // open file for output
 FILE *f = fopen("/tmp/out.txt","w");
 Candidate[BLOCK-1]=0;
 // repeats until the last tested prime exceed 10^9
 while (Candidate[BLOCK-1] < 1e9){
   //generate a block of primes
   for (uint32_t i=0; i<BLOCK; i++)
     Candidate[i] =  primesieve_next_prime(&it);
   // test them in parallel
#pragma omp parallel for schedule(dynamic)
   for (uint32_t i=0; i<BLOCK; i++)
     Res[i]=(Candidate[i]-2)%5 < 2 && Check(Candidate[i]);
   // save the results  
   for (uint32_t i=0; i<BLOCK; i++)
     if (Res[i] && Candidate[i]<1e9)
       fprintf(f,"%u\n",Candidate[i]);
 }
 fclose(f);
 return 0;
}

[Daniel Mondot]

Correction:

Out of 16000 terms of A064414, we have about:

- 37.8% that are primes of the form 20k+3

- 37.8% that are primes of the form 20k+7

- 24% that are composites of the form 24k+14

The sequence for the remaining 0.4% start: 1, 2, 4, 6, 9, 27, 49, 81, 98, 243, 343, 529, 729, 1849, 2187, 2401, 3698, 4489, 4802, 6561, 6889...etc

(note that the value of 0.4% is only true for 16000 terms, these numbers become increasingly rare)

Except for 1, 2 and 6, these are all powers of terms in A000057, or twice a power of a prime number (p^k or 2*p^k, where k>1, and p is a prime from A000057)

Also, the numbers of the form 2*p^k are numbers such that 2*p(k/2) ==14 mod 24

Note: because 2*2 is not == 14 mod 24, then 2*2^2=8 is not in A064414. and therefore no multiples of 8 are in the sequence.

But besides powers of 2 which stop at 4, I believe all powers of the rest of the primes in A000057 should be in A064414

Sorry about the mistake in my original email.

And do you think all this information and conjectures should be A064414?

Daniel.

[Sean A. Irvine]

> General questions: 
> Would it be worth having 5 new sequences:

> 1) The composite numbers of A064414 (I currently have about 3000 terms, hoping to get 10000 by the end of the week).

> 2) The first differences of A064414,

> 3) The first differences of its prime terms (A000057),

> 4) The first difference of its composite terms (new sequence (1)) ?

> 5) And since about 3/4 of the terms of A064414 are prime, and most of what remains are multiple of 24, we could have a sequence with just the remaining exceptions: the terms of A064414 that are neither primes or multiple of 24.

Yes. I see no problem adding all of these sequences.