１変数多項式の展開プログラム（いろんなサンプルコードつき）

expansion.rb

# １変数多項式の展開プログラム

class Polynomial 
	@@DEFAULT_MAX_ORDER = 200
	@@max_order = @@DEFAULT_MAX_ORDER   # 打ち切らずに計算させる最大の次数 ( (MAX_ORDER+1) 以上の次数は打ち切られる )

	attr_reader :coefficients, :variable


	# コンストラクタ
	def initialize coefficients_hash, variable="x"
		@coefficients = coefficients_hash
		@variable = variable
	end

	# 打ち切らずに計算させる最大の次数を設定する ( (MAX_ORDER+1) 以上の次数は打ち切られる )
	def self.setMaxOrder max_order
		@@max_order = max_order
	end


	# 他の多項式との積をとる
	def *(other)
		self.product(other)
	end


	# 他の多項式との積をとる
	def product other
		c = Hash.new

		a = @coefficients
		b = other.coefficients

		a_sorted = a.sort
		b_sorted = b.sort

		dim_c = a_sorted.max.first + b_sorted.max.first
		dim_c = (dim_c < @@max_order) ? dim_c : @@max_order

		(0..dim_c).each do |dim|
			# 次数がdimになる組み合わせ
			result = 0

			a.each_key do |k|
				l = dim - k

				if b.has_key?(l) then 
					result += a[k] * b[l]
				end
			end

			if result != 0 then
				c.store(dim, result)
			end

		end

		Polynomial.new c, @variable
	end


	# TeX 形式で出力
	def to_tex(sort=1)
		c = @coefficients
		if sort == 1 then
			# 昇順
			c = @coefficients.sort
		elsif sort == -1 then
			# 降順
			c = @coefficients.sort {|a, b| b <=> a }
		end

		str = ""

		c.each_with_index do |tuple, index|
			key   = tuple[0]
			value = tuple[1]

			if index > 0 then
				str += " "
			end

			if value > 0 then
				if index > 0 
					str += "+ "
				end
			elsif value < 0 then
				str += "- "
			end

			if value.abs == 1 then
				if key == 0 then					
					str += "1"
				elsif key != 1 then
					str += "#{@variable}^\{#{key}\}"
				else
					str += "#{@variable}"			
				end
			else
				if key == 0 then					
					str += "#{value.abs}"
				elsif key != 1 then
					str += "#{value.abs}#{@variable}^\{#{key}\}"
				else
					str += "#{value.abs}#{@variable}"					
				end
			end

		end

		str
	end

	# 係数配列として出力
	def to_array 
		array = Array.new

		c = @coefficients
		max_key = c.max[0]

		(max_key+1).times do |index|
			if c.has_key? index then
				array.push c[index]
			else
				array.push 0				
			end
		end

		array
	end

end



########################################################################
########################################################################
########################################################################
if($0==__FILE__) then
 	#テストコード

	require 'prime'

	###### <CASE 0>
	### シンプルな展開計算
	# (x^{2} + 2x + 1) (x^{2} + 4x + 4) = x^{4} + 6x^{3} + 13x^{2} + 12x + 4
	def testCase0
		a = Polynomial.new({2=>1, 1=>2, 0=>1})   # x^2 + 2x + 1
		b = Polynomial.new({2=>1, 1=>4, 0=>4})   # x^2 + 4x + 4
		c = a * b

		# 降べきの順に並べる
		puts "(#{a.to_tex(-1)}) (#{b.to_tex(-1)}) = #{c.to_tex(-1)}"

	end


	###### <CASE 1>
	### 二次形式 6xx+xy+yy に対応する保型形式 f のフーリエ係数
	def testCase1
		Polynomial.setMaxOrder(300)

		f = Polynomial.new({1=>1}, "q")
		(1..300).each do |n|
			if n%10 == 0 then
				puts "#{n}"
			end

			a1 = Polynomial.new({0=>1, n=>-1})
			a2 = Polynomial.new({0=>1, (23*n)=>-1})

			f = f * a1
			f = f * a2
		end

		puts ""
		puts f.to_tex
		puts ""
		p f.to_array
		puts ""


		### a(p) = 2 になる素数 p のリスト
		primes = Array.new
		f.to_array.each_with_index do |a_p, index|
			if index > 1 && Prime.prime?(index) && a_p == 2 then
				primes.push index
			end
		end
		p primes

	end


	###### <CASE 2>
	### 楕円曲線 y^2 = x^3 - x に対する保型形式 f のフーリエ係数
	def testCase2
		Polynomial.setMaxOrder(200)
		f = Polynomial.new({1=>1}, "q")                   # f := (1-q) 
		(1..200).each do |n|
			if n%10 == 0 then
				puts "#{n}"
			end
			a1 = Polynomial.new({0=>1, (4*n)=>-1}, "q")   # a1 := (1-q^{4n})
			a2 = Polynomial.new({0=>1, (8*n)=>-1}, "q")   # a2 := (1-q^{8n})
			f = f * a1*a1*a2*a2
		end
		puts ""
		puts f.to_tex
	end


	###### <CASE 3>
	### ラマヌジャンの Δ のフーリエ係数
	def testCase3
		Polynomial.setMaxOrder(100)
		delta = Polynomial.new({1=>1}, "q")          # delta := (1-q)
		(1..100).each do |n|
			if n%10 == 0 then
				puts "#{n}"
			end

			#puts "#{n}"
			a1 = Polynomial.new({0=>1, n=>-1}, "q")  # a1 := (1-q^n)
			24.times do 
				delta = delta * a1
			end
		end

		# 標準出力に掃き出してテスト
		puts ""
		puts delta.to_tex
		puts ""
		p delta.to_array


		# ラマヌジャン予想の確認
		File.open('ramanujan.dat', 'w') do |io|
			delta.to_array.each_with_index do |value, index|
				k = index
				if Prime.prime? k then
					io.puts "#{k} #{value.abs}"
				end
			end
		end

		puts ""

		##### Gnuplot の描画スクリプト
		puts "### try following gnuplot script"
		puts "$ gnuplot"
		puts "set logscale y"
		puts "plot \"ramanujan.dat\", 2*(x**(5.5))"

	end


	# テストケースの実行
	#testCase0
	#testCase1
	#testCase2
	testCase3

end ## if($0==__FILE__)