#!/usr/bin/env perl

# author: bliako
# date: 02/02/2019
# for:  https://perlmonks.org/?node_id=1229234
#
# Usage: mm.pl power [number-of-trials]
# Matrix is hard-coded.
#
# It calculates a matrix m to power p as efficiently as the author managed.
# The brute-force solution is : m x m x ... x m for p times.
# But, there are gains to harvest
# if we re-use results from intermediate multiplications.
# For example, m^14 = m^12 x m^2
# So, first calculate m^2 = m x m and cache it (that's 1 multiplication).
# Then calculate m^4 = m^2 x m^2 and cache it (that's another one).
# Then calculate m^6 = m^4 x m^2 and cache it (thats' another one).
# Then calculate m^8 = m^4 x m^4 and cache it (thats' another one).
# Then calculate m^12 = m^8 x m^4 and cache it (thats' another one).
# Then calculate m^14 = m^12 x m^2 and cache it (thats' another one).
# ...
# Then multiply the cached results with 1 more multiplications.
# Done with 7 multiplications!
# (as opposed to p-1=13 required by the brute-force approach)
# Now, there is scope for even more optimisation using
# the decomposition m^14 = m^4 x m^4 x m^4 x m^2
# which calculates the m^14 using only 5 multiplications!
# So, there is a naive randomised algorithm which picks plans
# at random and choses the one with the fewer number of multiplications
# (without doing the actual multiplications)
# That's the "plan" which is then returned as a code-ref which is
# executed to yield the matrix to the asked power using PDL's matrix mult (x)
#
# NOTE: m^a x m^b = m^(a+b)
#
use strict;
use warnings;

use PDL;
use PDL::Matrix;
use List::Util qw/shuffle/;

my $p = $ARGV[0]; # power to raise
my $n = $ARGV[1] // 500; # number of trials
# the matrix
my $m = PDL::Matrix->pdl([[3,5],[2,4]]);
if( $p > 65 ){
	print STDERR "$0 : power is too large, so I cowardly refuse to do any matrix multiplications, only the plan will be reported.\n";
	make_plan($p, $n);
} else {
	my $res = make_plan($p, $n)->($m);
	print "$0 : ". $m . " to the power of $p is".$res;
}
print "$0 : done.\n";

# Given a power (P), it tries different sets
# of even numbers whose sum equals P:
# e.g. P=14 yields 2,2,2,2,2,2,2 or 4,2,6,2 or ...
# Odd-P has just a 1 at the end.
# Then it tries each decomposition for how many multiplications would
# require and returns a code-ref to execute the plan given an input PDL
# matrix.
# No matrix multiplications below, just trying different decompositions
# with a naive randomised algorithm (it shuffles the items and picks the 1st 2)
# There is definetely scope for improvement, possibly with a genetic algorithm?
sub	make_plan {
	my $P = $_[0];
	my $trials = $_[1] // 500;

	my $best_plan = undef;
	my $bestTstr = undef;
	my $min_num_multiplications = 1000000;
	my ($sum, $nT);
	while( --$trials > 0 ){
		my @T = (2) x int($P/2);
		my $Tstr = join(",", @T);
		my $num_multiplications = 1; # the m^2 case

		my %cache = ();
		my $aplan = '$cache{2} = $m x $m;'."\n";
		while(1){
			@T = split(/,/, $Tstr);
			if( ($nT=scalar(@T)) == 1 ) { last }
			my @shuffledT = List::Util::shuffle(@T);
			my ($T0, $T1) = ($shuffledT[0], $shuffledT[1]);
			$sum = $T0 + $T1;
			$Tstr =~ s/${T0},${T1}/$sum/g;
			$Tstr =~ s/,+/,/g; $Tstr =~ s/^,//; $Tstr =~ s/,$//;
			if( ! exists $cache{$sum} ){
				$aplan .= '$cache{'.$sum.'} = $cache{'.$T1.'} x $cache{'.$T0.'};'."\n";
				$cache{$sum} = 1;
				$num_multiplications++;
			}
		}
		$num_multiplications += scalar(@T)-1;
		if( $num_multiplications < $min_num_multiplications ){
			$best_plan = $aplan; $bestTstr = $Tstr;
			print "make_plan() : found better plan with $num_multiplications multiplications: $Tstr\n";
			$min_num_multiplications = $num_multiplications;
		}
		print "make_plan() : trial $trials found $num_multiplications multiplications, best is $min_num_multiplications ...\n" if $trials % 10 == 0;
	}
	$best_plan .= 'my $result = ';
	$best_plan .= '$cache{'.$_.'} x ' foreach split /,/, $bestTstr;
	if( $P % 2 ){
		$min_num_multiplications++;
		$best_plan .= '$m;'."\n";
	} else { $best_plan =~ s/ x $/;\n/ }
	my @bestT = split /,/, $bestTstr;
	print "\nmake_plan() : final plan with $min_num_multiplications multiplications.\n";
	print "make_plan() : worst case requires ".($P%2?(int($P/2)+1):int($P/2))." multiplications and the naive case ".($P-1).".\n";
	print "make_plan() : the plan:\n$best_plan\n";

	$best_plan = 'my $m=$_[0]; my %cache=(); '.$best_plan.'return $result';
	my $ret = eval("sub{ $best_plan }");
	die "ERROR in eval of\nbegin code>>>>\n$best_plan\n<<<end code\n\n" if $@;
	return $ret
}
1;