#!/usr/bin/perl -w

# ----------------------------------------------------------------
# John Kerl
# kerl.john.r@gmail.com
# 2005-01-14
#
# This is a data generator for testing FFT implementations.  Example usage:
#
# bash$ prandgen 10
#      2.35684784579
#      1.11812729397
#      1.36974352416
#      0.59127890051
#      1.08356106794
#     -0.17648811137
#      0.25026955208
#      0.73697473761
#     -1.09526173095
#      0.51304069032
# ----------------------------------------------------------------

# Setup information:
# (1) Obtain my PMATLIB.pm;
# (2) Put PMATLIB.pm somewhere, e.g. the $HOME/bin directory;
# (3) Include that directory in the PERLLIB environment variable.
#     For bash, if PERLLIB exists:         export PERLLIB=$HOME/bin
#     For bash, if PERLLIB does not exist: export PERLLIB=$PERLLIB:$HOME/bin
#     For csh, if PERLLIB exists:          setenv PERLLIB $HOME/bin
#     For csh, if PERLLIB does not exist:  setenv PERLLIB ${PERLLIB}:$HOME/bin
use PMATLIB;

$do_seed     = 0;
$seed        = 0;
$do_complex  = 0;
$mean        = 0.0;
$stddev      = 1.0;
$num_rows    = 0;
$num_cols    = 0;
$symmetric   = 0;
$skew_symmetric = 0;
$pi          = 4.0 * atan2(1.0, 1.0);

# Save-variables for gauss_unit_distrib().
$ugd_iset = 0;
$ugd_gset = 0.0;

# xxx $do_uniform, $lo, $hi

sub usage {
	die
		"Usage: $0 [options] {num samples}\n" .
		"Or:    $0 [options] {num rows} {num cols}\n" .
		"Options:\n" .
		"  -r             Real output.\n" .
		"  -c             Complex output.\n" .
		"  -seed (seed)   Specified pseudorandom seed.\n" .
		"  -m    (mean)   Specified mean.  Default 0.0.\n" .
		"  -t    (stddev) Specified standard deviation.  Default 1.0.\n" .
		"  -symm          Produces a symmetric matrix.  Requires square dimensions.\n" .
		"  -sksymm        Produces a skew-symmetric matrix.  Requires square dimensions.\n";
	exit(1);
}

while (@ARGV > 1) {
	last unless $ARGV[0] =~ m/^-/;

	if ($ARGV[0] eq "-seed") {
		shift @ARGV;
		usage() unless @ARGV;
		$seed = shift @ARGV;
		$do_seed = 1;
	}
	elsif ($ARGV[0] eq "-m") {
		shift @ARGV;
		usage() unless @ARGV;
		$mean = shift @ARGV;
	}
	elsif ($ARGV[0] eq "-t") {
		shift @ARGV;
		usage() unless @ARGV;
		$stddev = shift @ARGV;
	}

	elsif ($ARGV[0] eq "-r") {
		$do_complex = 0;
		shift @ARGV;
	}
	elsif ($ARGV[0] eq "-c") {
		$do_complex = 1;
		shift @ARGV;
	}

	elsif ($ARGV[0] eq "-symm") {
		$symmetric      = 1;
		$skew_symmetric = 0;
		shift @ARGV;
	}

	elsif ($ARGV[0] eq "-sksymm") {
		$symmetric      = 0;
		$skew_symmetric = 1;
		shift @ARGV;
	}

	else {
		usage();
	}
}

if (@ARGV == 2) {
	$num_rows = shift @ARGV;
	$num_cols = shift @ARGV;
}
elsif (@ARGV == 1) {
	$num_rows = shift @ARGV;
	$num_cols = 1;
}
else {
	usage();
}

if ($symmetric && ($num_rows != $num_cols)) {
	usage();
}

if ($skew_symmetric && ($num_rows != $num_cols)) {
	usage();
}

if ($do_seed) {
	srand($seed);
}

if ($do_complex) {
	for ($i = 0; $i < $num_rows; $i++) {
		for ($j = 0; $j < $num_cols; $j++) {
			$mag = gauss_distrib($mean, $stddev);
			$phz = uniform_distrib(-$pi, $pi);
			$re = $mag * cos($phz);
			$im = $mag * sin($phz);
			printf " %18.11f %18.11f", $re, $im;
		}
		print "\n";
	}
}
else {
	if ($symmetric) {
		for ($i = 0; $i < $num_rows; $i++) {
			for ($j = 0; $j <= $i; $j++) {
				$A[$i][$j] = gauss_distrib($mean, $stddev);
				$A[$j][$i] = $A[$i][$j];
			}
		}
	}
	elsif ($skew_symmetric) {
		for ($i = 0; $i < $num_rows; $i++) {
			$A[$i][$i] = 0.0;
		}
		for ($i = 0; $i < $num_rows; $i++) {
			for ($j = 0; $j < $i; $j++) {
				$A[$i][$j] = gauss_distrib($mean, $stddev);
				$A[$j][$i] = -$A[$i][$j];
			}
		}
	}
	else {
		for ($i = 0; $i < $num_rows; $i++) {
			for ($j = 0; $j < $num_cols; $j++) {
				$A[$i][$j] = gauss_distrib($mean, $stddev);
			}
		}
	}
	print_matrix(\@A, $num_rows, $num_cols);
}

# ----------------------------------------------------------------
sub gauss_distrib {
	my ($mean, $stddev) = @_;
	return $mean + $stddev * gauss_unit_distrib();
}

# ----------------------------------------------------------------
# From _Numerical Recipes in C_.
# Returns a normally distributed deviate with zero mean and
# unit variance, using rand() as the source of uniform deviates.

sub gauss_unit_distrib {
	my ($fac, $rsq, $v1, $v2);
	if ($ugd_iset == 0) {
		# We don't have an extra deviate handy, so pick two uniform
		# numbers in the unit square, repeating until we get a pair
		# in the unit circle.
		do {
			$v1 = 2.0 * rand() - 1.0;
			$v2 = 2.0 * rand() - 1.0;
			$rsq = $v1*$v1 + $v2*$v2;
		} while ($rsq >= 1.0 || $rsq == 0.0);
		$fac = sqrt(-2.0 * log($rsq) / $rsq);
		# Now make the Box-Muller transformation to get two normal deviates.
		# Return one and save the other for next time.
		$ugd_gset = $v1 * $fac;
		$ugd_iset = 1;
		return $v2 * $fac;
	}
	else {
		$ugd_iset = 0;
		return $ugd_gset;
	}
}

# ----------------------------------------------------------------
sub uniform_distrib {
	my ($lo, $hi) = @_;
	return $lo + ($hi - $lo) * rand();
}