#!/usr/local/bin/perl
# saisho2jo.pl
#   n次の最小２乗法で近似曲線の係数、誤差を出す
#                  Jan 21 1999 by Tomonori NOMOTO

$usage = "Usage:  saisho2jo.pl ([options] [number])
Options: -d [fitする関数の次数] 次数を設定します
         -i [y切片の数値]       y切片(x=0時のyの値)を設定します";
unless(defined @ARGV){print $usage;}
$dim = 3;			# n次の最小2乗法
foreach $opt (@ARGV){
    if(defined $para){
	if($para == 1 && $opt =~ /^\d+/){$dim = $opt;}
	elsif($para == 2 && $opt =~ /^\d+/){$mode = 1;$intercept = $opt;}
	else{
	    print $usage;
	    exit;}
	undef $para;
	next;
    }
    if($opt =~ /^\-d/){$para=1;next}
    elsif($opt =~ /^\-i/){$para=2;next}
    $file = $opt;
}

open(DATA,$file);
$num = 0;				#データ数nの初期値
while(<DATA>){
    ($x[$num],$y[$num]) = split;
    $num++;
}
close DATA;

$r=$dim*2;			# 累乗数 x**? y**?
$rz=$dim;			# 累乗数 xy**?
for($a=0;$a<$num;$a++){
    $x = $x[$a];
    for($b=0;$b<=$r;$b++){
	$sx[$b] = $sx[$b] + $x ** $b;
    }
    $y = $y[$a];
    for($b=0;$b<=$rz;$b++){
	$syx[$b] = $syx[$b] + $y * $x ** $b;
    }
}

# Δを求める
for($r=0;$r<=$dim;$r++){
    $z = 1+$dim*2-$r;
    for($l=0;$l<=$dim;$l++){
	$z = $z - 1;
	$matrix[$l][$r]=$sx[$z];
    }
}
if($mode == 1){
    $delta=&determinant($dim - 1); 
}else{
    $delta=&determinant($dim);
}


for($c=0;$c<=$dim;$c++){
    for($r=0;$r<=$dim;$r++){
	if($r==$c){
	    $b = 1 + $dim;
	    for($l=0;$l<=$dim;$l++){
		$b = $b - 1;
		if($mode == 1){
		    $matrix[$l][$r] = $syx[$b] - $intercept * $sx[$b];
		}else{
		    $matrix[$l][$r]=$syx[$b];
		}

	    }
	}else{
	    $b = 1 + $dim * 2  - $r;
	    for($l=0;$l<=$dim;$l++){
		$b = $b - 1;
		$matrix[$l][$r]=$sx[$b];
	    }
	}
    }
    if($mode == 1){
	$const[$c]=&determinant($dim -1) / $delta;
    }else{
	$const[$c]=&determinant($dim) / $delta;
	 }
}
if($mode == 1){
    $const[$dim] = $intercept;
}

# 不偏分散u**2を求める
$syf=0;
for($a=0;$a<$num;$a++){
    $c=$dim+1;
    $fx=0;
    for($b=0;$b<=$dim;$b++){
	$c=$c-1;
	$fx=$fx+$const[$c]*$x[$a]**$b;
    }
    $syf=$syf +($y[$a]-$fx)**2;
}
$u = $syf/($num-$dim);

# 各係数の誤差
for($c=0;$c<=$dim;$c++){
    for($r=0;$r<=$dim;$r++){
	if($r==$c){
	    $b = 1 + $dim;
	    for($l=0;$l<=$dim;$l++){
		if($c==0){
		    $sig=1;
		}else{
		    $matrix[$l][$c]=$matrix[$l][0];
		    $sig=-1;
		}
		$b = $b - 1;
		$matrix[$l][0]=$syx[$b];
	    }
	}else{
	    $b = 1 + $dim * 2  - $r;
	    for($l=0;$l<=$dim;$l++){
		$b = $b - 1;
		$matrix[$l][$r]=$sx[$b];
	    }
	}
    }
    $dtwice=&error($dim,$sig);
    $err2[$c]=$dtwice*$u/($delta**2);
}

print " データ数n: $num\n";
print "      u**2= $u\n";
for($z=0;$z<=$dim;$z++){
    $mm = $dim - $z;
    print "x**$mmの係数: $const[$z]\n";
    $err[$z]=sqrt($err2[$z]);
    print "x**$mmの誤差: $err[$z]\n";
}

# 誤差計算の偏微分の２乗
sub error{
    undef @coefa;
    undef @coefb;
    $twice=0;
    ($d,$sig1)=@_;
    for($x=0;$x<=$d;$x++){for($y=0;$y<=$d;$y++){
	$backup[$x][$y]=$matrix[$x][$y];
    }}
    for($h=0;$h<=$d;$h++){		# 行
    for($x=0;$x<=$d;$x++){for($y=0;$y<=$d;$y++){
	$matrix[$x][$y]=$backup[$x][$y];
    }}
	if($h==$d){$sig2=1;}else{$sig2=-1;}
	    for($i=0;$i<=$d;$i++){ # 使う行を最後に持っていく
		$temp=$matrix[$d][$i];
		$matrix[$d][$i]=$matrix[$h][$i];
		$matrix[$h][$i]=$temp;
	    }
	for($i=0;$i<=$d;$i++){ # 先頭列を最後に持っていく
	    $temp=$matrix[$i][$d];
	    $matrix[$i][$d]=$matrix[$i][0];
	    $matrix[$i][0]=$temp;
	}
	$s=$d-$h;
	$coefa[$s] = $sig1*$sig2*&determinant($d-1);
    }
    for($h=0;$h<=$d;$h++){
	for($i=0;$i<=$d;$i++){
	    $j=$h+$i;
	    $coefb[$j]=$coefb[$j]+$coefa[$h]*$coefa[$i];
	}
    }
    for($h=0;$h<=$d*2;$h++){
	$twice=$twice+$coefb[$h]*$sx[$h];
    }
    return $twice;
}
# 行列式を求める
sub determinant{
    $nn = @_[0]; # $nはn+1行n+1列の行列のn
    &modify_determinant($nn);
    $det = 1;
    for($a=0;$a<=$nn;$a++){
	$det = $det * $matrix[$a][$a];
    }
    return $det;
}

# 行列式を変形する
sub modify_determinant{ 
    $n = @_[0]; # $nはn+1行n+1列の行列のn
    $msig=1;
    if($n == 0){
    }else{
	for($a=$n;$a>=0;$a--){
	    if($a == $n){
		$max = abs($matrix[$a][$n]);
		$maxl = $a;
		for($dc=$n - 1;$dc>=0;$dc--){
		    $mabs = abs($matrix[$dc][$n]);
		    if($mabs >= $max){
			$max = $matrix[$dc][$n];
			$maxl = $dc;
		    }
		}
		if($maxl != $a){
		    for($db=0;$db<=$a;$db++){
			$temp[$db] = $matrix[$maxl][$db];
			$matrix[$maxl][$db] = $matrix[$a][$db];
			$matrix[$a][$db] = $temp[$db];
		    }
		    $msig = $msig * -1;
		}
		$n_n = $matrix[$a][$n];
	    }else{
		$devide = $matrix[$a][$n] / $n_n;
		for($b=$n;$b>=0;$b--){
		    $matrix[$a][$b]=$matrix[$a][$b]-$matrix[$n][$b]*$devide;
		}
	    }
	}
	&modify_determinant($n - 1);
    }
}

# 行列を表示する (for Debug)
sub showmatrix{
    $matdim = @_[0];
    for($x=0;$x<=$matdim;$x++){
	for($y=0;$y<=$matdim;$y++){
	    print "$matrix[$x][$y]\t"
	    }
	print "\n";
    }
}