Devuan fork of gpsd
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

#### 120 lines 2.7 KiB Raw Permalink Blame History

 ``````#! /bin/perl # # From: Reuben Settergren # # I wrote a perl script to randomly generate matrices with smallish integers, # and test for a suitable determinant (typically +/-1), and ouput the matrix # and its inverse. # # There are switches to make it -s[ymmetric] and/or -d[iagonally strong], # which are necessary and good properties for covariance matrices. You may # need to grab Math::MatrixReal and Math::Factor::XS from cpan. use Math::MatrixReal; use Math::Factor::XS qw(prime_factors); use Getopt::Std; \$opt{n} = 1; getopts('hsdn:', \%opt); \$usage .= "easyinverse.pl [-hsd] [-n N]\n"; \$usage .= " -h this message\n"; \$usage .= " -s symmetric\n"; \$usage .= " -d diagonally dominant\n"; \$usage .= " -n N number to generate (DEF 1)\n"; if (\$opt{h}) { print \$usage; exit; } sub randint{ # but not 0 my \$min = shift or -5; my \$max = shift or +5; my \$ansa = 0; while (\$ansa == 0) { my \$r = rand; \$r *= (\$max - \$min + 1); my \$i = int(\$r); \$ansa = \$min + \$i; } return \$ansa; } \$size = 4; sub printmateq { my \$lbl = shift; my \$mat = shift; my \$round = shift; if (\$round) { for my \$r (1..\$size) { for my \$c (1..\$size) { my \$elt = \$mat->element(\$r,\$c); my \$pelt = abs(\$elt); my \$prnd = int(\$pelt * \$round + 0.5) / \$round; my \$rnd = \$prnd * (\$elt < 0 ? -1 : 1); \$mat->assign(\$r, \$c, \$rnd); }} } my @rows; for my \$r (1..\$size) { my @row = map {\$mat->element(\$r, \$_)} (1..\$size); my \$rowstr = '{' . (join ',', @row) . '}'; push @rows, \$rowstr; } my \$matstr = '{' . (join ',', @rows) . '}'; \$matstr =~ s!\.?000+\d+!!g; print "\$lbl = \$matstr;\n"; } for (1..\$opt{n}) { # requested number of matrices \$M = new Math::MatrixReal(\$size,\$size); while (1) { for \$r (1..\$size) { if (\$opt{d}) { \$M->assign(\$r,\$r, randint(5,10)) } else { \$M->assign(\$r,\$r, randint(-3, 3)) } for \$c (\$r+1..\$size) { \$i = randint(-5,5); \$M->assign(\$r, \$c, \$i); \$M->assign(\$c, \$r, (\$opt{s} ? \$i : randint(-5,5))); } } \$det = int(\$M->det + 0.5); next if \$det == 0; @factors = prime_factors(abs(\$det)); \$all_good_factors = 1; # cross your fingers! for \$f (@factors) { if (\$f != 2 || \$f != 5) { \$all_good_factors = 0; last; } } next unless \$all_good_factors; # test for positive semidefinite ~ all positive eigenvalues \$evals = \$M->sym_eigenvalues; \$all_positive_evals = 1; # keep hope alive! for \$r (1..\$size) { if (\$evals->element(\$r,1) <= 0) { \$all_positive_evals = 0; last; } } last if \$all_positive_evals; } printmateq(".mat", \$M); printmateq(".inv", \$M->inverse, 1000); } ``````