[ruby-list:21295] Re: class Polynomial

  豊福です。

  polynomial.rb の factorize を optimize しました。
パッチをメールの最後につけておきます。２割程度速く
なるようです。
  setPoly の方の修正はほとんどスピードには効いて
ないようですが、読みやすくはなったと思います。

  数式処理に関しては

    岩波講座 応用数学 方法９  計算代数と計算幾何

という本をみつけました。因数分解方法やグレブナ基底
の計算方法などものっていました。


  さて、もう一つの本題。
  実行の profile をとってみたのですが、Array#each
が２割くらいしめています。ですが polynomial.rb の
コードには each はほとんどなくて、調べてみると
each_with_index とfor i in 0..degD という部分が
Array#each の回数に加算されているようです。
  each_with_index の方は理解できますが、
for i in 0..degD の方は Range#each ですよね。
簡単なテストプログラムで for in 0..degD の profile
をとってみるとこちらは Range#each と出るところを
見るとバグでしょうか。
  ruby のバージョンは 1.4.3 です。
---
			豊福
			[e-mail]

*** polynomial.rb.orig	Tue Feb  1 03:33:53 2000
--- polynomial.rb	Fri Mar 10 16:03:15 2000
***************
*** 530,550 ****
  
  
  def  checkDivZp (dividend,divisor,prime) # true if divisible in Zp
! r=dividend.clone; degR=r.degree # Remainder(dividend)
! degD=divisor.degree; topDm=prime-Number.inv(divisor.array[degD],prime) # Divisor
! while degR>=degD
! 	q1=(r.array[degR]*topDm)%prime;
! 	dq=degR-degD
! 	if q1 != 0;
  		for i in (0..degD-1)
! 			r.array[i+dq] += divisor.array[i]*q1
  		end
  	end
- 	r.array[degR]=0
  	degR -= 1
  end
! r.normalize!
! r.array.each{|x| if (x%prime != 0) ;return false;end}
  return true
  end
  
--- 530,567 ----
  
  
  def  checkDivZp (dividend,divisor,prime) # true if divisible in Zp
! degR=dividend.degree
! degD=divisor.degree
! if degR < degD
! 	return false
! end
! darray=divisor.array
! rarray=dividend.array.clone # Remainder(dividend)
! topDm=prime-Number.inv(darray[degD],prime) # Divisor
! while degR>degD
! 	q1=(rarray[degR]*topDm)%prime;
! 	if q1!=0
! 		dq=degR-degD
  		for i in (0..degD-1)
! 			rarray[i+dq] += darray[i]*q1
  		end
  	end
  	degR -= 1
  end
! q1=(rarray[degR]*topDm)%prime;
! if q1!=0
! 	for i in (0..degD-1)
! 		if ((rarray[i] + darray[i]*q1) % prime) != 0
! 			return false
! 		end
! 	end
! else
! 	for i in (0..degD-1)
! 		if (rarray[i] % prime) != 0
! 			return false
! 		end
! 	end
! end
  return true
  end
  
***************
*** 671,683 ****
  
  
  
! def checkZp(i,prime)
  if i == 0;return false;end
  flg=false
- degD=DpolyM.degree
  while true
- 	if PolyN.degree<degD*2;return true;end
- 	if ! Number.checkDivZ?(PolyN.array[0],i,prime);return flg;end
  	if ! checkDivZp(PolyM,DpolyM,prime);return flg;end
  	if not resume(degD,prime);return flg;else flg=true;end
  	return flg # no need for loop now
--- 688,697 ----
  
  
  
! def checkZp(i,prime,degD)
  if i == 0;return false;end
  flg=false
  while true
  	if ! checkDivZp(PolyM,DpolyM,prime);return flg;end
  	if not resume(degD,prime);return flg;else flg=true;end
  	return flg # no need for loop now
***************
*** 685,714 ****
  end
  
  
! 
! def setPoly(d,top,prime)
  i=0; c=(prime/2).to_i
  while true
  	DpolyM.array[d]=i
  	if (d<=0);
- 		if top;return;end
  		if (i>0)&&Number.checkDivZ?(PolyN.array[0],i,prime)&&
! 			checkZp(i,prime)&&(PolyN.degree<DpolyM.degree*2);
  			throw(:setPolyTag)
  		end
! 	else if top&&(i>=1)
! 			if (PolyN.degree<DpolyM.degree*2);throw(:setPolyTag);end
! 			if Number.checkDivZ?(PolyN.array[PolyN.degree],i,prime)
! 				setPoly(d-1,false,prime)
! 			end
! 		else
! 			setPoly(d-1,top,prime)
  		end
  	end
- 	#i=i+1; if i>=prime;return;end
- 	if top; i=i+1; if i>=prime;return;end
- 	else if i<=c;i=prime-i-1; if i==c;return;end; else i=prime-i; end
- 	end
  end
  end
  
--- 699,729 ----
  end
  
  
! def setPoly2(d,prime,degD)
  i=0; c=(prime/2).to_i
  while true
  	DpolyM.array[d]=i
  	if (d<=0);
  		if (i>0)&&Number.checkDivZ?(PolyN.array[0],i,prime)&&
! 			checkZp(i,prime,degD)&&(PolyN.degree<degD*2);
  			throw(:setPolyTag)
  		end
! 	else 
! 		setPoly2(d-1,prime,degD)
! 	end
! 	if i<=c;i=prime-i-1; if i==c;return;end; else i=prime-i; end
! end
! end
! 
! def setPoly(d,prime)
! for degD in 1..d
! 	for i in 1..prime-1
! 		DpolyM.array[degD]=i
! 		if (PolyN.degree<degD*2);throw(:setPolyTag);end
! 		if Number.checkDivZ?(PolyN.array[PolyN.degree],i,prime)
! 			setPoly2(degD-1,prime,degD)
  		end
  	end
  end
  end
  
***************
*** 761,768 ****
  #printf "Set prime = %d\n",prime
  
  PolyM.array.replace(PolyN.modulo(prime).array);
! DpolyM.array.fill(0,0..degN)
! catch(:setPolyTag) do; setPoly(degN/2,true,prime);end
  if (PolyN.degree>0); Factor.push(PolyN);end
  
  # 最高次を正に揃える. お節介？
--- 776,783 ----
  #printf "Set prime = %d\n",prime
  
  PolyM.array.replace(PolyN.modulo(prime).array);
! DpolyM.array.clear
! catch(:setPolyTag) do; setPoly(degN/2,prime);end
  if (PolyN.degree>0); Factor.push(PolyN);end
  
  # 最高次を正に揃える. お節介？
***************
*** 773,779 ****
  
  module_function :factorize
  module_function :checkZp,:checkDivZp
! module_function :setPoly,:resume,:setVal,:setPolyR
  module_function :setV0Tbl,:setV1Tbl,:checkV
  
  end # module  Factorization
--- 788,794 ----
  
  module_function :factorize
  module_function :checkZp,:checkDivZp
! module_function :setPoly,:setPoly2,:resume,:setVal,:setPolyR
  module_function :setV0Tbl,:setV1Tbl,:checkV
  
  end # module  Factorization