spec/math/sieve_spec.cr

View this file on GitHub
Last update: 2022-01-02 17:13:44+09:00

Depends on

Code

require "spec"
require "../../src/math/sieve"
require "../../src/math/sieve_fast"

private N = 10000
private P = (2..N).select { |x| (2..Math.sqrt(x)).none? { |f| x % f == 0 } }

private def check_prime_division(x : Int32, division : Array({Int32, Int32}))
  division.all? { |factor, count| P.includes?(factor) && count >= 1 } &&
    division.each_cons(2, reuse: true).all? { |(t1, t2)| t1[0] < t2[0] } &&
    division.reduce(1) { |acc, (factor, count)| acc * factor**count } == x
end

{% for klass in [Sieve, SieveFast] %}
  describe {{klass}} do
    sieve = {{klass}}.new N
    it "#prime?" do
      (1..20).each do |x|
        sieve.prime?(x).should eq P.includes?(x)
      end
      expect_raises(ArgumentError) { sieve.prime? 0 }
      expect_raises(ArgumentError) { sieve.prime? -1 }
      expect_raises(ArgumentError) { sieve.prime? N + 1 }
    end

    it "#primes" do
      sieve.primes.should eq P
    end

    it "#prime_division" do
      (1..N).each do |x|
        result = sieve.prime_division(x).to_a
        check_prime_division(x, result).should be_true
      end
      expect_raises(ArgumentError) { sieve.prime_division(0) }
      expect_raises(ArgumentError) { sieve.prime_division(-1) }
      expect_raises(ArgumentError) { sieve.prime_division(N + 1) }
    end

    it "#each_factor" do
      (1..N).each do |x|
        result = [] of {Int32, Int32}
        sieve.each_factor(x) do |factor, count|
          result << {factor, count}
        end
        check_prime_division(x, result).should be_true
      end
      expect_raises(ArgumentError) { sieve.each_factor(0) { } }
      expect_raises(ArgumentError) { sieve.each_factor(-1) { } }
      expect_raises(ArgumentError) { sieve.each_factor(N + 1) { } }
    end

    it "#number_of_divisors" do
      (1..N).each do |x|
        expect = (1..x).count { |i| x % i == 0 }
        sieve.number_of_divisors(x).should eq expect
      end
      expect_raises(ArgumentError) { sieve.number_of_divisors(0) }
      expect_raises(ArgumentError) { sieve.number_of_divisors(-1) }
      expect_raises(ArgumentError) { sieve.number_of_divisors(N + 1) }
    end

    it "#sum_of_divisors" do
      (1..N).each do |x|
        expect = (1..x).select { |i| x % i == 0 }.sum(0i64)
        sieve.sum_of_divisors(x).should eq expect
      end
      expect_raises(ArgumentError) { sieve.sum_of_divisors(0) }
      expect_raises(ArgumentError) { sieve.sum_of_divisors(-1) }
      expect_raises(ArgumentError) { sieve.sum_of_divisors(N + 1) }
    end
  end
{% end %}

require "spec"

# require "../../src/math/sieve"
class Sieve
  getter size : Int32, factor : Array(Int32), primes : Array(Int32)

  def initialize(@size)
    @factor = Array(Int32).new(@size + 1, 0)
    @primes = [] of Int32
    sqrt_size = Math.sqrt(@size).to_i + 1
    (2..size).each do |x|
      next unless @factor[x] == 0
      @factor[x] = x
      @primes << x
      next if sqrt_size < x
      (x * x).step(to: size, by: x) do |y|
        @factor[y] = x if @factor[y] == 0
      end
    end
  end

  def prime?(x : Int) : Bool
    raise ArgumentError.new unless 1 <= x <= size
    factor[x] == x
  end

  private class PrimeDivisionIterator
    include Iterator({Int32, Int32})

    def initialize(@current : Int32, @factor : Array(Int32))
    end

    def next
      return stop if @current == 1
      element = @factor[@current]
      count = 0
      while @current != 1 && @factor[@current] == element
        count += 1
        @current //= element
      end
      {element, count}
    end
  end

  def prime_division(x : Int)
    raise ArgumentError.new unless 1 <= x <= size
    PrimeDivisionIterator.new(x, factor)
  end

  def each_factor(x : Int, &) : Nil
    raise ArgumentError.new unless 1 <= x <= size
    while x > 1
      element = @factor[x]
      count = 0
      while x != 1 && @factor[x] == element
        count += 1
        x //= element
      end
      yield(element, count)
    end
  end

  def number_of_divisors(x : Int) : Int32
    raise ArgumentError.new unless 1 <= x <= size
    cnt = 1
    each_factor(x) do |_, c|
      cnt *= c.succ
    end
    cnt
  end

  def sum_of_divisors(x : Int) : Int64
    raise ArgumentError.new unless 1 <= x <= size
    sum = 1i64
    each_factor(x) do |elem, cnt|
      sum *= (elem.to_i64 ** cnt.succ - 1) // elem.pred
    end
    sum
  end
end

# require "../../src/math/sieve_fast"
class SieveFast
  getter size : Int32, factor : Array(Int32), primes : Array(Int32)

  private def value(i)
    2 * i + 3
  end

  private def index(x)
    (x - 3) // 2
  end

  # index :  0  1  2  3  4  5  6  7  8  9 10 11
  # value :  3  5  7  9 11 13 15 17 19 21 23 25
  # factor:  3  5  7  3 11 13  3 17 19  3 23  5

  private def sift(n)
    # factor(i)       = value(i)          = 2i + 3
    # index_square(i) = index(value(i)^2) = 2i^2 + 6i + 3
    # factor(i + 1)       - factor(i)       = 2
    # index_square(i + 1) - index_square(i) = 4i + 8 = factor(i) + factor(i + 1)
    i, factor, index_square = 0, 3, 3
    while index_square < n
      if @factor[i] == 0
        @factor[i] = factor
        @primes << factor
        index_square.step(to: n - 1, by: factor) do |j|
          @factor[j] = factor if @factor[j] == 0
        end
      end
      i += 1
      index_square += factor
      factor += 2
      index_square += factor
    end
    while i < n
      if @factor[i] == 0
        @factor[i] = factor
        @primes << factor
      end
      i += 1
      factor += 2
    end
  end

  def initialize(@size)
    @primes = [2]
    n = (@size - 1) // 2
    @factor = Array(Int32).new(n, 0)
    sift(n)
  end

  def prime?(x : Int32)
    raise ArgumentError.new unless 1 <= x <= size
    if x.even?
      x == 2
    elsif x == 1
      false
    else
      @factor[index(x)] == x
    end
  end

  private class PrimeDivisionIterator
    include Iterator({Int32, Int32})

    @done2 : Bool

    def initialize(@current : Int32, @factor : Array(Int32))
      @done2 = @current.odd?
    end

    def next
      return stop if @current == 1

      unless @done2
        @current //= 2
        count = 1
        while @current.even?
          @current //= 2
          count += 1
        end
        @done2 = true
        return {2, count}
      end

      element = @factor[(@current - 3) // 2]
      @current //= element
      count = 1
      while @current != 1 && @factor[(@current - 3) // 2] == element
        count += 1
        @current //= element
      end
      {element, count}
    end
  end

  def prime_division(x : Int32)
    raise ArgumentError.new unless 1 <= x <= size
    PrimeDivisionIterator.new(x, @factor)
  end

  def each_factor(x : Int32, & : Int32, Int32 ->) : Nil
    raise ArgumentError.new unless 1 <= x <= size
    if x.even?
      count = 1
      x //= 2
      while x.even?
        count += 1
        x //= 2
      end
      yield 2, count
    end
    while x > 1
      element = @factor[index(x)]
      x //= element
      count = 1
      while x != 1 && @factor[index(x)] == element
        count += 1
        x //= element
      end
      yield element, count
    end
  end

  def number_of_divisors(x : Int) : Int32
    cnt = 1
    each_factor(x) do |_, c|
      cnt *= c + 1
    end
    cnt
  end

  def sum_of_divisors(x : Int) : Int64
    sum = 1i64
    each_factor(x) do |elem, cnt|
      sum *= (elem.to_i64 ** cnt.succ - 1) // (elem - 1)
    end
    sum
  end
end

private N = 10000
private P = (2..N).select { |x| (2..Math.sqrt(x)).none? { |f| x % f == 0 } }

private def check_prime_division(x : Int32, division : Array({Int32, Int32}))
  division.all? { |factor, count| P.includes?(factor) && count >= 1 } &&
    division.each_cons(2, reuse: true).all? { |(t1, t2)| t1[0] < t2[0] } &&
    division.reduce(1) { |acc, (factor, count)| acc * factor**count } == x
end

{% for klass in [Sieve, SieveFast] %}
  describe {{klass}} do
    sieve = {{klass}}.new N
    it "#prime?" do
      (1..20).each do |x|
        sieve.prime?(x).should eq P.includes?(x)
      end
      expect_raises(ArgumentError) { sieve.prime? 0 }
      expect_raises(ArgumentError) { sieve.prime? -1 }
      expect_raises(ArgumentError) { sieve.prime? N + 1 }
    end

    it "#primes" do
      sieve.primes.should eq P
    end

    it "#prime_division" do
      (1..N).each do |x|
        result = sieve.prime_division(x).to_a
        check_prime_division(x, result).should be_true
      end
      expect_raises(ArgumentError) { sieve.prime_division(0) }
      expect_raises(ArgumentError) { sieve.prime_division(-1) }
      expect_raises(ArgumentError) { sieve.prime_division(N + 1) }
    end

    it "#each_factor" do
      (1..N).each do |x|
        result = [] of {Int32, Int32}
        sieve.each_factor(x) do |factor, count|
          result << {factor, count}
        end
        check_prime_division(x, result).should be_true
      end
      expect_raises(ArgumentError) { sieve.each_factor(0) { } }
      expect_raises(ArgumentError) { sieve.each_factor(-1) { } }
      expect_raises(ArgumentError) { sieve.each_factor(N + 1) { } }
    end

    it "#number_of_divisors" do
      (1..N).each do |x|
        expect = (1..x).count { |i| x % i == 0 }
        sieve.number_of_divisors(x).should eq expect
      end
      expect_raises(ArgumentError) { sieve.number_of_divisors(0) }
      expect_raises(ArgumentError) { sieve.number_of_divisors(-1) }
      expect_raises(ArgumentError) { sieve.number_of_divisors(N + 1) }
    end

    it "#sum_of_divisors" do
      (1..N).each do |x|
        expect = (1..x).select { |i| x % i == 0 }.sum(0i64)
        sieve.sum_of_divisors(x).should eq expect
      end
      expect_raises(ArgumentError) { sieve.sum_of_divisors(0) }
      expect_raises(ArgumentError) { sieve.sum_of_divisors(-1) }
      expect_raises(ArgumentError) { sieve.sum_of_divisors(N + 1) }
    end
  end
{% end %}