page top

package k2;

;# -----------------------------------------------------------------------------
;# K2/LIB::PERL LIBRARY++                                                      .
;#                                                                             .
;#       Mathematics Library                                                   .
;#                                                                             .
;#                             All Right Reserved, CopyRight (C) 2008 by Dr.K  .
;# -----------------------------------------------------------------------------


;#//////////////////////////////////////////////////////////////////////////////
;# Perl Package Module Import
;#//////////////////////////////////////////////////////////////////////////////
 #use strict;
;#//////////////////////////////////////////////////////////////////////////////


;#//////////////////////////////////////////////////////////////////////////////
;# Prototype
;#//////////////////////////////////////////////////////////////////////////////
;# 乱数関係の関数
;# BinomRand($p, $n)
;# NormalRand($m, $s, $c)

# 素数関係の関数 x 4
;# HashCode($s)
;# isPrime($v)
;# PrimeNear($v, $f)
;# PrimeFactor($v, $f)

;# 計算用関数 x 14
;# Log2($v, $e)
;# Power($v, $p)
;# Total(@dat)
;# Av(@dat)
;# Vx(@dat)
;# Sd(@dat)
;# Se(@dat)
;# t0(@dat)
;# Sd1(@dat)
;# Variance($sw, @dat)
;# Permutation($n, $r)
;# Combination($n, $r)
;# Homogeneous($n, $r)
;# Factorial($n);

;# Radian対応三角関数 x 10
;# Pai($n)
;# Rad2Deg($r)
;# Deg2Rad($d)
;# Sin($d)
;# Cos($d)
;# Tan($d)
;# Cot($d)
;# Asin($d)
;# Acos($d)
;# Atan2($x, $y)

;# 数値変換 x 6
;# cnv10To2($v)
;# cnv2To10($v)
;# cnv10To16($v)
;# cnv16To10($v)
;# rgb($col)
;# color($R, $G, $G)

;# Total = 34
;#//////////////////////////////////////////////////////////////////////////////


;#//////////////////////////////////////////////////////////////////////////////
;# 乱数関係の関数
;#//////////////////////////////////////////////////////////////////////////////


sub PaiRnd{
	my ($n) = @_;
	my $r = 0;
	for(1..$n){$r++ if(rand() ** 2 + rand() ** 2 <= 1)}
	return($r);
}

# 二項分布乱数
# $p=成功確率
# $n=試行回数
sub BinomRand{
	my ($p, $n) = @_;

	my $r = 0;
	# ベルヌーイ試行
	for(1..$n){$r++ if($p > rand())}

    return($r);
}

# 正規分布乱数
# $m=平均値(mean)
# $s=標準偏差(sigma)
# $c=カウント(count)
sub NormalRand{
	my ($m, $s, $c) = @_;

	my ($r1, $r2, $rr);

	$r1 = rand();
	$r2 = rand();

	if($c % 2){
		$rr = $s * sqrt(-2 * log($r1)) * sin(2 * 3.14159 * $r2) + $m;
	}else{
		$rr = $s * sqrt(-2 * log($r1)) * cos(2 * 3.14159 * $r2) + $m;
	}

	return($rr);
}


;#//////////////////////////////////////////////////////////////////////////////
;# 素数関係の関数
;#//////////////////////////////////////////////////////////////////////////////


# ハッシュ関数(素数法)
sub HashCode{
	my ($s) = @_[0];

	my $h = 0;
	my $ret = '';
	my $ls = length($s);
	for(my $i = 0; $i < $ls; $i++){
		$h +=  37 * ord(substr($s, $i, 1));
		#$ret .= $h % 1987;
	}
	$ret =  PrimeFactor($h, 1);
	return($ret);
}

# 素数の判定(高速)
sub isPrime{
	my ($v) = @_;
	return 0 if ($v < 2);
	return 1 if ($v == 2);
	return 0 if ($v % 2 == 0);
	for(my $i = 3; $i * $i <= $v; $i += 2){
		return(0) unless($v % $i);
	}
	return(1);
}

# 近隣素数の取得
# $f=1::-側、$f=2::+側
sub PrimeNear{
	my ($v, $f) = @_;
	my $expr = 0;
	$f = MinMax($f, 1, 2);
	$f = ($f == 1) ? -1 : 1;
	for(my $i = $v - 1; $i > 0; $i = $i + $f){
		if(isPrime($i)){
			$expr = $i;
			last;
		}
	}
	return($expr);
}

# 素因数分解
# $f=1::素因数、$f=2::素数(低速)
# 引数に '+, -, *, **' の演算子が使用出来る
sub PrimeFactor{
	my ($v, $f) = @_;
	$f = 1 unless($f);
	$f = MinMax($f, 1, 2);

	# use only operator of '+, -, *, **'
	# 式の展開($v = 1000 || $v = 10,000 || $v = 10 * 2 + 1)
	$v = eval($v) if($v =~ /[1-9]\d* *(\+|-|\*|\*\*) *[1-9]\d*/);

	# 値の制限($v > 0 && $v < 999999999)
	return('') if($v !~ /^[1-9]\d{0,8}$/);

	my ($ret, %result);
	my @prime = (2, 3);

	for(my $i = 0; $prime[$i] ** 2 <= $v; $i++){
		unless($v % $prime[$i]){
			$result{$prime[$i]}++;
			$v /= $prime[$i];
			$i--;
			next;
		}
		if($i == $#prime){
			my $p = $prime[-1] + 2;
			for(my $j = 1; $prime[$j] ** 2 <= $p; $j++){
				unless($p % $prime[$j]){
					$p += 2;
					$j = 1;
					redo;
				}
			}
			push(@prime, $p);
		}
	}

	if($f == 1){
		if(%result){
			$result{$v}++;
			my $expr = join ' x ',
					map{$result{$_} == 1 ? "$_" : "$_ ^ $result{$_}"}
					sort{$a <=> $b} keys %result;
			$ret = $expr;
		}else{
		}
	}else{
		$ret = (%result) ? 0 : 1;
	}

	return($ret);

	# 正規表現素因数分解
	#my ($ret, $n);
	#for($n = ('o' x $v); $n =~ /^(oo+?)\1+$/; $n =~ s/$1/o/g){
	#    $ret .= length($n).'('.length($1).')'.',';
	#}
}


;#//////////////////////////////////////////////////////////////////////////////
;# 計算用関数
;#//////////////////////////////////////////////////////////////////////////////


# 底'$e'が指定できる自然対数
sub Log2{
	my ($v, $e) = @_;
	if($e){
		# 底が'$e'
		return(log($v) / log($e));
	}else{
		# 底はネイピア数(e)
		return(log($v));
	}
}

# トータル計算
sub Total{
	my (@dat) = @_;
	my $tl = 0;
	foreach(@dat){$tl += $_}
	return($tl);
}

# 平均値計算
# Mean
sub Av{
	my (@dat) = @_;
	return(0)unless(@dat);
	my $av = Total(@dat);
	$av /= @dat;
	return($av);
}

# 分散値計算
# Variance
sub Vx{
	my (@dat) = @_;
	return(Variance(1, @dat));
}

# 標準偏差計算
# Standard Deviation
sub Sd{
	my (@dat) = @_;
	return(Variance(2, @dat));
}

# 標準誤差計算
# Standard Error
sub Se{
	my (@dat) = @_;
	return(0)unless(@dat);
	return(Sqr(Vx(@dat) / @dat));
}

# 検定統計量(t0)計算
# chisq test
sub t0{
	my ($d1, $d2) = @_;
	return(0)unless($d1.$d2);
	my $t0 = abs(Av(@$d1) - Av(@$d2)) / Sqr(Se(@$d1) + Se(@$d2));
	return($t0);
}

# ($sw=1)分散
# ($sw=2)標準偏差           _
# 不偏分散(s^2 = 1/2��(xi - x)^2)
sub Variance{
	my ($sw, @dat) = @_;

	my $ret = 0;
	my $num = @dat;
	$sw = 1 unless($sw);
	if($num){
		my $ave = Av(@dat);

		while(@dat){
			my $v = shift(@dat);
			$ret += Power(($v - $ave), 2);
		}
		$ret = $ret / $num;
		if($sw == 2){$ret = Sqr($ret)}
	}
	return($ret);
}

# 標準偏差1(違う方法::少し答えが違う::大きい)
sub Sd1{
	my (@dat) = @_;
	my $sum1 = 0;
	my $sum2 = 0;

	foreach my $x(@dat){
		$sum1 += $x;
		$sum2 += Power($x, 2);
	}

	my $nm = @dat;
	my $sd = sqrt( ($sum2 - Power($sum1, 2) / $nm) / ($nm - 1) );
	return ($sd);
}

# 順列::nPr
sub Permutation{
	my ($n, $r) = @_;
	my $p = 1;
	for(my $i = 1; $i <= $r; $i++){
		$p *= ($n - $i + 1);
	}
	return($p);
}

# 組合わせ::nCr
# recursive call(再帰呼び出し)
# 整数演算に除算を使うは好ましくない
sub Combination{
	my ($n, $r) = @_;

	# Combination1
	if($n == $r){
		return(1);
	}elsif($r == 1){
		return($n);
	}elsif($n < $r){
		return(0);
	}else{
		return(Combination($n - 1, $r) + Combination($n - 1, $r - 1));
	}

	# Combination2
	# 上の関数の方がかなり早い
	# for(my $i = 1; $i <= $r; $i++){
	#    $c = ($c * ($n - $i + 1) / $i);
	#};
	# $c = Round($c);
}

# 重複組合わせ::nHr
sub Homogeneous{
	my ($n, $r) = @_;
	my $h = Combination($r + $n - 1, $r);
	return($h);
}

# 階乗::n!
sub Factorial{
	my ($n) = @_;
	my $p = 1;
	for(my $i = 1; $i <= $n; $i++){
		$p *= $n;
	}
	return($p);
}


;#//////////////////////////////////////////////////////////////////////////////
;# 三角関数
;#//////////////////////////////////////////////////////////////////////////////


# 円周率の定数宣言
use constant MTH_PAI_1 => 3.1415926535;
use constant MTH_PAI_2 => MTH_PAI_1 * 2;
#use constant DEG_1 => MTH_PAI_2 / 360;

# πを得る
# Pai()->π、Pai(2)->2π
sub Pai{
	my ($n) = @_;
	$n = 1 unless($n);
	return(MTH_PAI_1 * $n);
}

# RadianからDegree変換
sub Rad2Deg{
	my ($r) = @_;
	$r = 360 * ($r / MTH_PAI_2);
	return($r);
}

# DegreeからRadian変換
sub Deg2Rad{
	my ($d) = @_;
	$d = $d * (MTH_PAI_2 / 360);
	return($d);
}

# Degree対応sin
sub Sin{
	my ($d) = @_;
	my $r = Deg2Rad($d);
	$r = sin($r);
	return($r);
}

# Degree対応cos
sub Cos{
	my ($d) = @_;
	my $r = Deg2Rad($d);
	$r = cos($r);
	return($r);
}

# Degree対応tan
sub Tan{
	my ($d) = @_;
	my $r = Deg2Rad($d);
	$r = sin($r) / cos($r);
	return($r);
}

# Degree対応ctan
sub Cot{
	my ($d) = @_;
	my $r = Deg2Rad($d);
	$r = ($r == 0) ? '' : 1 / tan($r);
	return($r);
}

# Degree対応asin
sub Asin{
	my ($d) = @_;
	my $r = Deg2Rad($d);
	$r = atan2($r, sqrt(1 - $r * $r));
	return($r);
}

# Degree対応acos
sub Acos{
	my ($d) = @_;
	my $r = Deg2Rad($d);
	$r = atan2(sqrt(1 - $r * $r), $r);
	return($r);
}

# Degree対応atan
sub Atan2{
	my ($x, $y) = @_;
	my $r = atan2($y, $x);
	$r = Rad2Deg($r);
	return($r);
}


;#//////////////////////////////////////////////////////////////////////////////
;# 数値変換
;#//////////////////////////////////////////////////////////////////////////////


# 10進->2進変換(Bin)
sub cnv10To2{
	my ($v) = @_;
	$v = sprintf("%b", $v);
	return($v);
}

# 2進->10進変換(Decimal)
sub cnv2To10{
	my ($v) = @_;
	my $v = oct("0b".$v);
	return($v);
}

# 10進->16進変換(Hex)
sub cnv10To16{
	my ($v) = @_;
	$v = sprintf("%x", $v);
	return($v);
}

# 16進->10進変換(Decimal)
sub cnv16To10{
	my ($v) = @_;
	$v = hex($v);
	return($v);
}

# ColorCodeをRGBに分解
# ([#]123456) -> (R, G, B)
sub rgb{
	my($col) = @_;
	my(@RGB, $i, $s);
	$col = '000000'.$col;
	$col = substr($col, length($col) - 6, 6);

	for($i = 0; $i < 6; $i += 2){
		$s = substr($col, $i, 2);
		$RGB[@RGB] = hex($s);
	}
	return(@RGB);
}

# RGBをColorCodeに合成
# (R, G, B) -> ([#]123456)
sub color{
	my ($R, $G, $B) = @_;
	my $col = '';
	foreach($R, $G, $B){
		$_ = unpack("H2", pack("C", $_));
		$col .= "$_";
	}
	return($col);
}

## K2/LIBPERL LIBRARY++ ------------------------------------------------------
1;                                                      ## Presented by Dr.k2 ##
## EOF -------------------------------------------------------------------------