Segmented Sieve Memory Weirdness

output.txt

n # => 100
31: 9.408 kiB
38: 9.696 kiB
44: 9.76 kiB
50: 9.84 kiB
54: 9.84 kiB
58: 9.84 kiB
62: 9.872 kiB
64: 9.872 kiB
62: 9.872 kiB
64: 9.872 kiB
62: 9.872 kiB
64: 9.872 kiB
62: 9.872 kiB
64: 9.872 kiB
62: 9.872 kiB
64: 9.872 kiB
58: 9.872 kiB
62: 9.904 kiB
64: 9.904 kiB
62: 9.904 kiB
64: 9.904 kiB
62: 9.904 kiB
64: 9.904 kiB
68: 9.904 kiB
76: 9.968 kiB
44: 9.984 kiB
50: 10.064 kiB
54: 10.064 kiB
58: 10.064 kiB
62: 10.096 kiB
64: 10.096 kiB
62: 10.096 kiB
64: 10.096 kiB
62: 10.096 kiB
64: 10.096 kiB
62: 10.096 kiB
64: 10.096 kiB
62: 10.096 kiB
64: 10.096 kiB
58: 10.096 kiB
62: 10.128 kiB
64: 10.128 kiB
62: 10.128 kiB
64: 10.128 kiB
62: 10.128 kiB
64: 10.128 kiB
62: 10.128 kiB
64: 10.128 kiB
58: 10.128 kiB
62: 10.16 kiB
64: 10.16 kiB
62: 10.16 kiB
64: 10.16 kiB
68: 10.16 kiB
76: 10.16 kiB
44: 10.176 kiB
50: 10.256 kiB
54: 10.256 kiB
58: 10.256 kiB
62: 10.288 kiB
64: 10.288 kiB
62: 10.288 kiB
64: 10.288 kiB
62: 10.288 kiB
64: 10.288 kiB
62: 10.288 kiB
64: 10.288 kiB
62: 10.288 kiB
64: 10.288 kiB
58: 10.288 kiB
62: 10.32 kiB
64: 10.32 kiB
62: 10.32 kiB
64: 10.32 kiB
62: 10.32 kiB
64: 10.32 kiB
58: 10.32 kiB
62: 10.352 kiB
64: 10.352 kiB
62: 10.352 kiB
64: 10.352 kiB
68: 10.352 kiB
76: 10.352 kiB
44: 10.368 kiB
50: 10.448 kiB
54: 10.448 kiB
58: 10.448 kiB
62: 10.48 kiB
64: 10.48 kiB
62: 10.48 kiB
64: 10.48 kiB
62: 10.48 kiB
64: 10.48 kiB
62: 10.48 kiB
64: 10.48 kiB
62: 10.48 kiB
64: 10.48 kiB
58: 10.48 kiB
62: 10.512 kiB
64: 10.512 kiB
62: 10.512 kiB
64: 10.512 kiB
62: 10.512 kiB
64: 10.512 kiB
58: 10.512 kiB
62: 10.544 kiB
64: 10.544 kiB
62: 10.544 kiB
64: 10.544 kiB
58: 10.544 kiB
62: 10.576 kiB
64: 10.576 kiB
62: 10.576 kiB
64: 10.576 kiB
68: 10.576 kiB
76: 10.688 kiB
44: 10.704 kiB
50: 10.784 kiB
54: 10.784 kiB
58: 10.784 kiB
62: 10.816 kiB
64: 10.816 kiB
62: 10.816 kiB
64: 10.816 kiB
62: 10.816 kiB
64: 10.816 kiB
62: 10.816 kiB
64: 10.816 kiB
62: 10.816 kiB
64: 10.816 kiB
58: 10.816 kiB
62: 10.848 kiB
64: 10.848 kiB
62: 10.848 kiB
64: 10.848 kiB
62: 10.848 kiB
64: 10.848 kiB
62: 10.848 kiB
64: 10.848 kiB
58: 10.848 kiB
62: 10.88 kiB
64: 10.88 kiB
62: 10.88 kiB
64: 10.88 kiB
58: 10.88 kiB
62: 10.912 kiB
64: 10.912 kiB
68: 10.912 kiB
76: 10.912 kiB
44: 10.928 kiB
50: 11.008 kiB
54: 11.008 kiB
58: 11.008 kiB
62: 11.04 kiB
64: 11.04 kiB
62: 11.04 kiB
64: 11.04 kiB
62: 11.04 kiB
64: 11.04 kiB
62: 11.04 kiB
64: 11.04 kiB
62: 11.04 kiB
64: 11.04 kiB
58: 11.04 kiB
62: 11.072 kiB
64: 11.072 kiB
62: 11.072 kiB
64: 11.072 kiB
62: 11.072 kiB
64: 11.072 kiB
58: 11.072 kiB
62: 11.104 kiB
64: 11.104 kiB
62: 11.104 kiB
64: 11.104 kiB
58: 11.104 kiB
62: 11.136 kiB
64: 11.136 kiB
62: 11.136 kiB
64: 11.136 kiB
68: 11.136 kiB
76: 11.136 kiB
44: 11.152 kiB
50: 11.232 kiB
54: 11.232 kiB
58: 11.232 kiB
62: 11.264 kiB
64: 11.264 kiB
62: 11.264 kiB
64: 11.264 kiB
62: 11.264 kiB
64: 11.264 kiB
62: 11.264 kiB
64: 11.264 kiB
62: 11.264 kiB
64: 11.264 kiB
58: 11.264 kiB
62: 11.296 kiB
64: 11.296 kiB
62: 11.296 kiB
64: 11.296 kiB
62: 11.296 kiB
64: 11.296 kiB
58: 11.296 kiB
62: 11.328 kiB
64: 11.328 kiB
62: 11.328 kiB
64: 11.328 kiB
58: 11.328 kiB
62: 11.36 kiB
64: 11.36 kiB
68: 11.36 kiB
76: 11.36 kiB
44: 11.376 kiB
50: 11.456 kiB
54: 11.456 kiB
58: 11.456 kiB
62: 11.488 kiB
64: 11.488 kiB
62: 11.488 kiB
64: 11.488 kiB
62: 11.488 kiB
64: 11.488 kiB
62: 11.488 kiB
64: 11.488 kiB
62: 11.488 kiB
64: 11.488 kiB
58: 11.488 kiB
62: 11.52 kiB
64: 11.52 kiB
62: 11.52 kiB
64: 11.52 kiB
62: 11.52 kiB
64: 11.52 kiB
62: 11.52 kiB
64: 11.52 kiB
58: 11.52 kiB
62: 11.552 kiB
64: 11.552 kiB
62: 11.552 kiB
64: 11.552 kiB
58: 11.552 kiB
62: 11.584 kiB
64: 11.584 kiB
68: 11.584 kiB
76: 11.584 kiB
44: 11.6 kiB
50: 11.68 kiB
54: 11.68 kiB
58: 11.68 kiB
62: 11.712 kiB
64: 11.712 kiB
62: 11.712 kiB
64: 11.712 kiB
62: 11.712 kiB
64: 11.712 kiB
62: 11.712 kiB
64: 11.712 kiB
62: 11.712 kiB
64: 11.712 kiB
58: 11.712 kiB
62: 11.744 kiB
64: 11.744 kiB
62: 11.744 kiB
64: 11.744 kiB
62: 11.744 kiB
64: 11.744 kiB
58: 11.744 kiB
62: 11.776 kiB
64: 11.776 kiB
62: 11.776 kiB
64: 11.776 kiB
58: 11.776 kiB
62: 11.808 kiB
64: 11.808 kiB
62: 11.808 kiB
64: 11.808 kiB
68: 11.808 kiB
76: 12.016 kiB
GC.stats.total_bytes.humanize_bytes # => "11.7kiB"

sieve_of_eratosthenes.cr

require "bit_array"

struct BitArray
  def fill(value : Bool) : self
    @bits.to_slice(malloc_size).fill (value ? ~0u32 : 0u32)
    self
  end
end

def log_total_bytes(line)
  print line, ": ", (GC.stats.total_bytes / 1000), " kiB\n"
end

def sieve_of_eratosthenes(upto n)
  possible_primes = BitArray.new n + 1, initial: true
  possible_primes[0] = false
  possible_primes[1] = false

  (2..).take_while { |i| i <= Math.sqrt(n) }.each do |i|
    if possible_primes[i]?
      (i**2).upto(n).step(i).each do |j|
        possible_primes[j] = false
      end
    end
  end

  possible_primes.each_index.select { |index| possible_primes[index] }.to_a
end

def segmented_sieve_of_eratosthenes(upto n)
  log_total_bytes __LINE__

  delta = Math.sqrt(n).to_i

  # known_primes = Array(Int32).new initial_capacity: (n / Math.log(n)).to_i
  known_primes = sieve_of_eratosthenes upto: delta

  log_total_bytes __LINE__

  # we know that, after the first instance, 2 will never be counted as prime, so we use it as a filler at the end
  # TODO: make sure this doesn't cause problems for small n
  possible_primes = BitArray.new delta
  (delta + 1).upto(n).step(delta).each do |segment_min|
    log_total_bytes __LINE__

    segment_max = segment_min + delta - 1
    prime_divisor_limit = Math.sqrt(segment_max)
    relevant_prime_divisors = known_primes.each.take_while { |prime| prime <= prime_divisor_limit }

    log_total_bytes __LINE__

    possible_primes.fill true

    log_total_bytes __LINE__

    pre_segment_min = segment_min - 1
    relevant_prime_divisors.each do |p|
      log_total_bytes __LINE__

      starting_index = p - (pre_segment_min % p) - 1
      starting_index.upto(delta - 1).step(p).each do |j|
        log_total_bytes __LINE__
        possible_primes[j] = false
        log_total_bytes __LINE__
      end
    end

    log_total_bytes __LINE__

    possible_primes.each_index do |index|
      if possible_primes[index]
        known_primes << (segment_min + index)
      end
    end

    log_total_bytes __LINE__
  end

  known_primes
end

if PROGRAM_NAME.ends_with? "sieve_of_eratosthenes"
  n = 100
  p! n
  sieve = segmented_sieve_of_eratosthenes(upto: n)
  p! GC.stats.total_bytes.humanize_bytes
end