On the Infinity of Twin Primes and other K-tuples

I just released my paper - On the Infinity of Twin Primes and other K-tuples - and thought it would be of interest here how I used Ruby|Crystal to compute the data presented in the paper.

https://www.academia.edu/41024027/On…other_K-tuples

https://www.scribd.com/document/4364…other-K-tuples

I used Ruby version 2.7.0-preview2 (final due Christmas 2019, as usual) to write my Ruby scripts. It has the nice Enumerable#tally method, which made life really simple for me. I also used my Rubygem primes-utils to generate the primes (which I have a working Crystal version of, but haven’t released yet).

So after installing primes-utis I can do oneliners as below to find all the prime gaps within any range p to p*p, as long as p isn’t too large.

require 'primes/utils'

2.7.0-preview1 :037 > n = 5; n.primes(n*n).each_cons(2).map{ |a,b| b-a }.tally.sort
 => [[2, 3], [4, 3]]

2.7.0-preview1 :038 > n = 7; n.primes(n*n).each_cons(2).map{ |a,b| b-a }.tally.sort
 => [[2, 4], [4, 5], [6, 2]]

2.7.0-preview1 :039 > n = 11; n.primes(n*n).each_cons(2).map{ |a,b| b-a }.tally.sort
 => [[2, 8], [4, 9], [6, 7], [8, 1]]

2.7.0-preview1 :040 > n = 13; n.primes(n*n).each_cons(2).map{ |a,b| b-a }.tally.sort
 => [[2, 9], [4, 11], [6, 10], [8, 1], [10, 1], [14, 1]]

2.7.0-preview1 :041 > n = 17; n.primes(n*n).each_cons(2).map{ |a,b| b-a }.tally.sort
 => [[2, 16], [4, 14], [6, 17], [8, 1], [10, 3], [12, 2], [14, 1]]

2.7.0-preview1 :042 > n = 19; n.primes(n*n).each_cons(2).map{ |a,b| b-a }.tally.sort
 => [[2, 17], [4, 17], [6, 19], [8, 1], [10, 5], [12, 2], [14, 3]]

2.7.0-preview1 :043 > n = 23; n.primes(n*n).each_cons(2).map{ |a,b| b-a }.tally.sort
 => [[2, 21], [4, 23], [6, 25], [8, 7], [10, 7], [12, 4], [14, 3]]

2.7.0-preview1 :044 > n = 29; n.primes(n*n).each_cons(2).map{ |a,b| b-a }.tally.sort
 => [[2, 29], [4, 30], [6, 38], [8, 12], [10, 15], [12, 7], [14, 4], [18, 1]] 

But to process the gaps for very large p, I had to get more creative to get past large memory usage, and time. So I first created a Ruby version to count residue gaps where I directly computed the residues in order, and then computed|counted the gap size between all consecutive residues. When I got it to work correctly, I translated it to Crystal for speed. Here’s the Crystal version.

require "big"

def gaps(p)
  h = Hash(UInt64, UInt64).new(0)
  if p == 2; h[2] = 1; return h end
  primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
  modpg = primes[0..primes.index(p)].product(1u64)
  rescnt = primes[0..primes.index(p)].map{ |p| p-1 }.product(1u64)
  puts "modp#{p} = #{modpg}"
  puts "rescnt = #{rescnt}"
  pc, inc = 5u64, 2
  #r0 = second = [primes.index(p)+1..primes.index(p)+1]
  #first = primes[primes.index(p)+2..primes.index(p)+2]
  #second, first = primes[primes.index(p)+1..primes.index(p)+2]
  #r0 = second
  first = 43u64
  second = r0 = 41u64
  limit = modpg >> 1
  while pc < limit
    if pc.gcd(modpg) == 1
       first = pc
       gap = first - second
       h[gap] += 1
       second = first
       print "\r #{(first * 100)/limit}%"
    end
    pc += inc; inc = inc ^ 0b110
  end
  puts
  h.delete(0)
  h = h.to_a.map{ |k,v| (k == 2 || k == 4) ? [k,2*v+1] : [k,2*v]}.to_h
  if h[r0 - 1] == nil; h[r0 - 1] = 2 else h[r0 - 1] += 2 end
  h.to_a.sort_by { |key, value| key }
end

puts gaps 37

A problem with this version is the commented out section creates compiler errors because it squawks raising nil errors for the index values I’m trying to get. So I just hard code the values in so I can run it and get my results.

This routine is what I use to generate the data in Fig 3. in the paper.

I’m currently running this Crystal script to generate all the residue gaps for P37. I started it on a 2011 Lenovo laptop with an I5, 4 core cpu, and 6 GB of ram. But since this version is single threaded the extra cores don’t help me. My real work computer is a System76, I7 - 8 threads laptop, which I have to do all my other stuff with.

I began running it on Sunday, November 10, 2019 at 10:45 am. At time of writing (atow), Saturday, November 23, 4:35 pm, it’s now 60%+ done, clocking in at about 1% every 5+ hours, with estimated finish in 21/22 days, or by November 30, 2019

I can go faster by splitting the residues range into segments and doing them in parallel. I’d like to try to do this in Crystal to give me an incentive to learn its new parallel programming powers (PPP), and produce the data for the next larger generator P41, which has 40 times more residues than P37 {(41 - 1) for p = 41.}

If anyone wants to play with it please let me know.
Any suggestions to speed it up, fix the indexing issue, and make the code more efficient, are welcome.

4 Likes

you can do x.index(...).not_nil! if you are sure the thing you are trying to find the index for is there.

Also, is that snippet the entire program? I ask because maybe there’s a way to optimize it? I’d like to benchmark it or profile it to see if it can be improved. I already see some temporary arrays being created though I’m not sure they are the bottleneck.

Yes, the code is the whole program.

Using your tip on the index issue, this works.

require "big"

def gaps(p)
  h = Hash(UInt64, UInt64).new(0)
  if p == 2; h[2] = 1; return h end
  primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
  i = primes.index(p).not_nil!
  modpg  = primes[0..i].product(1u64)
  rescnt = primes[0..i].map{ |p| p-1 }.product(1u64)
  puts "modp#{p} = #{modpg}"
  puts "rescnt = #{rescnt}"
  pc, inc, r0 = 5u64, 2, 0
  second, first = primes[i+1..i+2]
  r0 = second
  limit = modpg >> 1
  while pc < limit
    if pc.gcd(modpg) == 1
       first = pc
       gap = first - second
       h[gap] += 1
       second = first
       print "\r #{(first * 100)/limit}%"
    end
    pc += inc; inc = inc ^ 0b110
  end
  puts
  h.delete(0)
  h = h.to_a.map{ |k,v| (k == 2 || k == 4) ? [k,2*v+1] : [k,2*v]}.to_h
  if h[r0 - 1] == nil; h[r0 - 1] = 2 else h[r0 - 1] += 2 end
  h.to_a.sort_by { |key, value| key }
end

puts gaps 23

And it compiles as: $ crystal build gaps.cr --release

And I run it as: $ time ./gaps
.
To make it more user friendly, I’d write it to take the value p as a cli input so I wouldn’t have to recompile for different p values

Great! I’ll play around with it later today if I can.

So, I profiled the code. Most of the time is spent in the printf line, also converting the progress number into a float.

I recommend either removing the progress report or showing it as an int but only when it changes:

  percent = 0
  while pc < limit
    if pc.gcd(modpg) == 1
      first = pc
      gap = first - second
      h[gap] += 1
      second = first
      new_percent = ((first * 100)/limit).to_i
      if new_percent != percent
        print "\r #{percent}%"
      end
      percent = new_percent
    end
    pc += inc; inc = inc ^ 0b110
  end

The times without progress report, or with that enhanced progress report, compared to the float report, vary from 6 seconds to 35 seconds.

So the algorithm seems to be fine! The progress report is the slow thing.

You can also use this to show elapsed time and remaining time:

  print "\rProgress: #{percent}%"
  time = Time.monotonic
  while pc < limit
    if pc.gcd(modpg) == 1
      first = pc
      gap = first - second
      h[gap] += 1
      second = first
      new_percent = (first * 100)/limit
      if (new_percent - percent) >= 0.001
        percent = new_percent
        elapsed = Time.monotonic - time
        remaining = (elapsed / new_percent) * (100 - new_percent)
        print "\rProgress: #{percent.round(3)}%, Elapsed: #{elapsed}, Remaining: #{remaining}"
      end
    end
    pc += inc; inc = inc ^ 0b110
  end

In my machine, for p 37 it seems it will take about 2 days to finish.

You can spawn a fiber that prints the progress each second:

spawn do
  loop do
    print "#{(first * 100)/limit}%\r"
    sleep 1
  end
end

Use the-Dpreview_mt, else nothing will be printed.

Please use Time.monotonic for measurements ;)

Below is optimized version, using the Prime Generator Sequence (pgs) for P7
(for those who read the paper) to more optimally generate the residues.

This code ran fine in Crystal 0.31.1 until I tried it for p = 31.

When compiled like: $ crystal build gaps.cr --release
it compiles, and runs with no apparent prolblem, until it finishes, then…

➜  Crystal Projects time ./gaps6b                     
modp31 = 200560490130
rescnt = 30656102400

Unhandled exception: Arithmetic overflow (OverflowError)
  from ???
  from gaps6b.cr:0:0 in 'gaps'
  from gaps6b.cr:43:6 in '__crystal_main'
  from /home/jzakiya/crystal/share/crystal/src/crystal/main.cr:47:14 in 'main'
  from __libc_start_main
  from ../sysdeps/x86_64/start.S:122:0 in '_start'
  from ???
./gaps6b  6537.84s user 0.07s system 99% cpu 1:48:58.21 total

I spent parts of 3 days doing everything I could think of to find what was
causing the overflow, to no avail. Finally, after remembering I could turn it
off at compile time, I recompiled like this.

$ $ crystal build -Ddisable_overflow gaps6b.cr --release

and it ran to completion and get these results

➜  Crystal Projects time ./gaps6b                                       
modp31 = 200560490130
rescnt = 30656102400
 98%
[{2, 1931585729}, {4, 1931585729}, {6, -311471742}, {8, -1770392096}, {10, -1309530646}, {12, -2088758088}, 
{14, 903350042}, {16, 423955224}, {18, 512670088}, {20, 106604280}, {22, 126682650}, {24, 88337252}, {26, 18298102}, 
{28, 16461600}, {30, 9169532}, {32, 833688}, {34, 1075458}, {36, 620632}, {38, 77042}, {40, 128988}, {42, 40636}, 
{44, 3516}, {46, 1794}, {48, 1296}, {50, 504}, {52, 20}, {54, 84}, {56, 12}, {58, 2}]
./gaps6b  6564.59s user 0.19s system 99% cpu 1:49:25.18 total

Looking at Fig 3. in my paper shows the first 6 results (for gaps 2-12) are incorrect, because those values are bigger than i32 bit values, but those variables are supposed to be u64. I really have no idea what going on or how to fix it. The correct values for the first six gaps should be:

{2, 6226553025}, {4, 6226553025}, {6, 8278462850},
{8, 2524575200}, {10, 2985436650}, {12, 2206209208},

Can someone point me to where the conversion to 64-bit values isn’t happening.

Also:

As a user, and to help new users get past this ANOMALY, I make these 2 suggestions.

  1. Make -Ddisable_overflow the refault compiler behavior,

  2. Explicitly state in the error message for overflow, to recompile with the disable flag.

I can’t tell you how annoying this was, and the fact there is no easy way for me to even
know where in the program the overflow is occurring. This will REALLY turn people off
from Crystal, especially for people doing number crunching with it, and coming from other languages.

require "big"

# $ crystal build -Ddisable_overflow gaps.cr --release

def gaps(p)
  h = Hash(UInt64, UInt64).new(0)
  if p == 2; h[2] = 1; return h end
  primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47] of UInt64
  pgspn = [2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4,  2, 4, 2, 
           4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10]
  pc = 11u64
  i = primes.index(p).not_nil!
  modpg  = primes[0..i].product
  rescnt = primes[0..i].map{ |p| p-1 }.product
  puts "modp#{p} = #{modpg}"
  puts "rescnt = #{rescnt}"
  second, first = primes[i+1..i+2]
  r0 = second
  percent = 0
  limit = modpg >> 1
  while pc < limit
    if pc.gcd(modpg) == 1
      first = pc
      gap = first - second
      h[gap] += 1
      second = first
      new_percent = ((first * 100)/limit).to_i
      if new_percent != percent
        print "\r #{percent}%"
      end
      percent = new_percent
    end
    inc = pgspn.shift
    pc += inc; pgspn << inc
  end
  puts
  h.delete(0)
  h = h.to_a.map{ |k,v| (k == 2 || k == 4) ? [k,2*v+1] : [k,2*v]}.to_h
  if h[r0 - 1] == nil; h[r0 - 1] = 2 else h[r0 - 1] += 2 end
  h.to_a.sort_by { |key, value| key }
end

puts gaps 31

EDITED:
From future discussion on the code producing overflow errors replace:

h = h.to_a.map{ |k,v| (k == 2 || k == 4) ? [k,2*v+1] : [k,2*v]}.to_h

with

h = h.to_a.map{ |k,v| (k == 2 || k == 4) ? [k,v*2+1] : [k,v*2]}.to_h

and compile as normal: $ crystal build gaps.cr --release

The type of an operation like a * b gets the type of the operator on the left.

I see you have 2*v+1. That means the type will be Int32, because we have 2 (Int32) multiplied by something else.

I’m almost sure that if you change that to be v*2 + 1 and v*2 then you will get a correct result, without overflow.

As a user, and to help new users get past this ANOMALY, I make these 2 suggestions.

As you already saw, you disabled overflow but got incorrect results. An overflow error means there’s a bug in your code. You definitely want to spot these and avoid them. So we will never go with either of the 2 options you propose.

In fact, it’s likely that in the future there won’t be any way to disable overflow, because that essentially means “I want potential bugs in my code”.

One thing we could do, however, is to change math operations so when you do a * b the type you get is the maximum type of both a and b, though that’s also a bit confusing, for example if you have a var x with type Int8 and you do x *= 2 then it won’t work anymore (because the result would be Int32 now, not Int8).

Another alternative is to restrict operations of a * b to exactly the same types. That means that x *= 2 in the above example will actually work thanks you automatic cast, but it won’t work if you do x *= y and y is of a different type.

So every choice has its pros and cons.

Can the compiler maybe spot possible overflow problems? and allow missmatched arithmetic operations if they are definitely safe, and complain if not?

Like multiply int16 and int8 and store in int32 or more, this will always be fine. but complain when target type has less bits.

No, I don’t see a clear way to implement something like that.

Thanks @asterite for the information.

However, this brings up a really important philosophical question about what you want Crystal to be, and for whom.

I’m coming primarily from Ruby, where Matz makes programmer happiness the number one priority. In Ruby, all this number casting isn’t the programmer’s problem it’s the VM. I don’t have to worry (in general) about what the types of numbers are, and certainly not the order they are used in for simple arithmetic.

Arithmetic should just work! The programmer shouldn’t have to act like a compiler and be forced to write code to satisfy it. The compiler should be written to make my life easy and do the obvious things.

The default understanding should be if a person writes a + b or x * y the intent is to get an answer that would make sense with a calculator, thus return the answer that fits into the largest number type. If the user wants something else, then they can explicitly cast the numbers to do so.

But at the very, very, least, this type of stuff needs to be repeatedly and vociferously documented and talked about. It’s not intuitive that I, as the programmer, when writing a Crystal program, need to know, and be always concerned with, simple commutative arithmetic operations like + and * not being type commutative. I don’t want to be forced to have that job.

So, the software compiler (writers) should write it to make us programmers happy.

Rust has realized this, and made its compiler give you nice suggestions to try when you, invariably, mess up (as you will, especially in the beginning) to not leave you totally frustrated and cursing it (which you’ll do anyway).

I’ve spent three (now 4) days trying to get this to work in Crystal (to speed up the Ruby version) because the compiler didn’t do the obvious thing and make 2*v return the max type value. This makes me sad and unhappy. :face_with_head_bandage:

We are aware of this problem (or at least I am aware of this problem). I’d like to improve the situation. The problem is… how? Going Ruby’s route, making a single numeric type and deciding at runtime what to do, converting things behind the hood, essentially means it will work like in Ruby: it will be intuitive and nice, but it will be slow. Having it like in Crystal means it will be fast but you have to be careful.

Do you have a suggestion of how to make it work well?

The default understanding should be if a person writes a + b or x * y the intent is to get an answer that would make sense with a calculator, thus return the answer that fits into the largest number type.

The only way to make that work would be to have a single number type that can hold arbitrary sizes, because if you use Int64 you can still overflow if you multiply two of them.

So again: what would be a solution we could use?

Finally, if your code is heavily numeric oriented I would recommend always using big_int (even in that Hash(UInt64, UInt64) thing) or using another language that’s more suited for math stuff, like Julia (though I don’t know whether they have overflow or not, or how they handle those things).

Look, I like Crystal, I don’t want to not use it, I want it to become more useful (to me).

One nice feature Ruby has now in its repl is a did_you_mean gem that gives Rust like suggestions when some common coding errors are made.

I’m not a compiler writer, so I don’t know the intricacies you have to contend with, but I know for some compilers/languages, if they see multiplication by 2 ,4, 8, etc they convert them to shift operators. I would think, if you have a constant at compile time that’s the easiest type to convert because it’s an explicit immutable value.

But even if you kept everything as is, compiler wise, to make programmers lives easier, publicize the quirks prominently in the docs, blogs, tutorials etc. For example, you could promote the rule, always have constant values be the right side operand, etc, and maintain that style in all the source code to reinforce it for users.

I can live with stuff I know about (and remember), it’s the unknown unknowns that kill programmers.

Oh by the way, changing 2*v to v*2 did correct the overflow error, compiled without the disable flag: $ crystal build gaps.cr --release

➜  Crystal Projects time ./gaps6b                    
modp31 = 200560490130
rescnt = 30656102400
 98%
[{2, 6226553025}, {4, 6226553025}, {6, 8278462850}, 
{8, 2524575200}, {10, 2985436650}, {12, 2206209208}, 
{14, 903350042}, {16, 423955224}, {18, 512670088}, 
{20, 106604280}, {22, 126682650}, {24, 88337252}, 
{26, 18298102}, {28, 16461600}, {30, 9169532}, 
{32, 833688}, {34, 1075458}, {36, 620632}, {38, 77042}, 
{40, 128988}, {42, 40636}, {44, 3516}, {46, 1794}, {48, 1296}, 
{50, 504}, {52, 20}, {54, 84}, {56, 12}, {58, 2}]
./gaps6b  6576.18s user 0.19s system 99% cpu 1:49:36.91 total

Cool!

I’m thinking, maybe the best way would be to have the compiler follow this rule:

Another alternative is to restrict operations of a * b to exactly the same types. That means that x *= 2 in the above example will actually work thanks you automatic cast, but it won’t work if you do x *= y and y is of a different type.

I think most of the annoyance before automatic type casting between numeric types (when a literal value is given) came from having to say things like 2_u64 when the compiler could obviously see that with 2 it was enough.

I might try to play with this idea and see how it affects the compiler’s source code and standard library. If it doesn’t bring that many issues (i.e.: if we need to make just a few modifications to make it work) then it might be worth discussing it in a GH issue to eventually push it forward.

That said, I still don’t know how we can handle the Int64 * Int64 possibly overflowing. The only way to solve that is to always use BigInt. But that’s not as friendly as in Ruby, I guess…

3 Likes

But @asterite the Int64 * Int64 case is (should) already be a concern for people doing numerical programming. That’s why in my case, where all the values are always positive, It was simpler/better to explicitly use UInt64.

When you know the domain size of your results, you can program to use the most efficient representation for them. When you don’t use the biggest possible, you need to customize your code to get correct results anyway.

There obviously can’t be a 100% perfect world (when dealing with infinity, or infinitesimals) but you can take care of 95% of the obvious and let programmers know what needs to be done to satisfy their special use cases. Sort of like Rails. It’s opinionated by default, but it’s all just Ruby underneath, so you can still do what you want.