2012年3月2日金曜日

桁の小さい平方根の計算

とりあえず桁の小さい平方根の計算の処理は完成しました
次はBIGFloatを利用してもう少し大きい桁の平方根を計算したいです
完成まで長かったwて、やってないから長いんですけどw

後、HSPでグラフを描くモジュールも作成したいです
Perlでグラフを描くのは僕の技量では無理w

use strict;
use warnings;

my ($r,$ans);

($r,$ans)=&small_sqrt(20000.0,0);
print $r." ".$ans."\n";

($r,$ans)=&small_sqrt(0.0002,0);
print $r." ".$ans."\n";

($r,$ans)=&small_sqrt(20000.0,1);
print $r." ".$ans."\n";

($r,$ans)=&small_sqrt(0.0002,1);
print $r." ".$ans."\n";

# 平方根の計算。計算する桁が小さい
sub small_sqrt{

# 引数A:sqrtしたい値
my $value=$_[0];
# 引数B:これが1だとprint debug
my $on_debug=$_[1];

# 返り値用
# $r:error情報。-1:数値が大きすぎる:-2数値が小さすぎる
# $ans:解
my ($r,$ans);

if ($value>=1.0){
# sqrtしたい値が1.0より大きい場合
($r,$ans)=&small_sqrt_L_in($value,$on_debug);
}else{
# sqrtしたい値が1.0より小さい場合
($r,$ans)=&small_sqrt_S_in($value,$on_debug);
}

return ($r,$ans);
}

# 平方根の計算。計算したい値は1.0より大きい
sub small_sqrt_L_in{

# 計算したい値
my $p=$_[0];
# これが1だとdebug
my $on_debug=0;
$on_debug=$_[1];

if ($on_debug==1){
printf "sqrt(".$p.")"."\n";
}

# &some_divはsqrtしたい値が1以上の時呼び出す
# sqrtしたい値は10の何乗の2乗か?求める
my ($r,$r_A,$r_B)=&some_div($p);

# &in_calc_0p1to0p0000001は実際に計算する
my $ans=&in_calc_0p1to0p0000001($r_A);

if ($on_debug==1){
printf $r."\n";
printf $r_A."\n";
printf $r_B."\n";

printf $ans*$r_B."\n";
}

return ($r,$ans*$r_B);

}

# 平方根の計算。計算したい値は1.0より小さい
sub small_sqrt_S_in{

# 計算したい値
my $p2=$_[0];
# これが1だとdebug
my $on_debug=0;
$on_debug=$_[1];

if ($on_debug==1){
printf "sqrt(".$p2.")"."\n";
}

# &some_mulはsqrtしたい値が1より小さい時呼び出す
# sqrtしたい値は0.1の何乗の2乗か?求める
my ($r2,$r_A2,$r_B2)= &some_mul($p2);

# &in_calc_0p1to0p0000001は実際に計算する
my $ans2=&in_calc_0p1to0p0000001($r_A2);

if ($on_debug==1){
printf $r2."\n";
printf $r_A2."\n";
printf $r_B2."\n";

printf $ans2*$r_B2."\n";
}

return ($r2,$ans2*$r_B2);

}

# sqrtしたい値が1以上の時呼び出す
# sqrtしたい値は10の何乗の2乗か?求める
sub some_div{
# $r:error通知用。普通は1。10の301以上だとerror
# $r_A:繰り返しで/100.0づつしていく
# $r_b:10のX乗の2乗か?のXの値を入れる
my ($r,$r_A,$r_B)=(1,$_[0],1.0);
# $r_A=&mul_10_X_count(301);
if ($r_A>&mul_10_300_count()){
$r=-1;
}
while ($r_A >= 100.0) {
$r_A/=100.0;
$r_B*=10.0;
}
return ($r,$r_A,$r_B);
}
# sqrtしたい値が1より小さい時呼び出す
# sqrtしたい値は0.1の何乗の2乗か?求める
sub some_mul{
# $r:error通知用。普通は1。10の301以上だとerror
# $r_A:繰り返しで*100.0づつしていく
# $r_b:0.1のX乗の2乗か?のXの値を入れる
my ($r,$r_A,$r_B)=(1,$_[0],1.0);
# $r_A=&mul_0p1_X_count(285);
if ($r_A<&mul_0p1_285_count()){
$r=-2;
}
while ($r_A <= 1.0) {
$r_A*=100.0;
$r_B/=10.0;
}
return ($r,$r_A,$r_B);
}
# 計算の処理。in_calcを8回呼び出す
sub in_calc_0p1to0p0000001{
# $p:この数値以下で最も大きい数値を導きたい。等しくてもいい
# $ret0~$ret8:計算結果が入る。少しづつ小さな単位を計算して代入
# &in_calcの3つめの引数:計算する単位
my $p=$_[0];
my $ret01=&in_calc($p,0,1);
my $ret02=&in_calc($p,$ret01,0.1);
my $ret03=&in_calc($p,$ret02,0.01);
my $ret04=&in_calc($p,$ret03,0.001);
my $ret05=&in_calc($p,$ret04,0.0001);
my $ret06=&in_calc($p,$ret05,0.00001);
my $ret07=&in_calc($p,$ret06,0.000001);
my $ret08=&in_calc($p,$ret07,0.0000001);
return $ret08;
}
# 計算の処理。一桁だけ計算
sub in_calc{
# この数値以下で最も大きい数値を導きたい。等しくてもいい
my $expression=$_[0];
# 計算結果を入れる
my $exp01=$_[1];
# 計算の単位
my $unit01=$_[2];
# 返り値用
my $ans01=0;
# 乗算の基
my @for_mul;
for (0..9){
$for_mul[$_]=$exp01+($_*$unit01);
}
# 乗算の基を2乗
my @of_mul;
for (0..9){
$of_mul[$_]=$for_mul[$_]*$for_mul[$_];
# print "i:$for_mul[$_]:i*i:$of_mul[$_]\n";
}
# X以下で大きい数値を求める
for (0..9){
if ($of_mul[$_]<=$expression){
$ans01=$for_mul[$_];
}
}
# print "解の候補:$ans01\n";
$ans01;
}
#printf &mul_10_X_count(0)."\n";
#printf &mul_10_X_count(3)."\n";
#printf &mul_10_300_count()."\n";
#printf &mul_10_300_count()*&mul_0p1_300_count()."\n";
#printf &mul_0p1_X_count(0)."\n";
#printf &mul_0p1_X_count(3)."\n";
#printf &mul_0p1_300_count()."\n";
#printf &mul_0p1_300_count()*&mul_10_300_count()."\n";

# 10の引数A乗の数を返す
sub mul_10_X_count{
my $x=1.0;
my $i=$_[0];
$i--;
for(0..$i){
$x*=10;
}
return $x;
}
# 10の300乗の数を返す
sub mul_10_300_count{
return &mul_10_X_count(300);
}

# 0.1の引数A乗の数を返す
sub mul_0p1_X_count{
my $x=1.0;
my $i=$_[0];
$i--;
for(0..$i){
$x/=10;
}
return $x;
}
# 0.1の285乗の数を返す
sub mul_0p1_285_count{
return &mul_0p1_X_count(285);
}