#!/usr/local/bin/perl
# -*- perl -*-
 
sub shortusage {
    print <<USAGE;
 
plotr - plot functions defined by recurrence equations
logplotr - same as plotr, but plots using logarithmic scaling
 
Usage:
  plotr -p -e -n nlimit -t tlimit fn fn ...
  logplotr -p -e -n nlimit -t tlimit fn fn ...
 
Options:
  -e    Print plot as Encapsulated Postscript to plotr.eps if normal exit.
  -h    Help.  Prints a longer explanation of plotr.
  -n    Specifies limit on value of N (x-axis).  Default is 100.
        More than 1000 can be very slow for certain types of recurrences.
  -p    Print plot as Postscript to plotr.ps if normal exit.
  -t    Specifies limit on T(N) (y-axis).  Default is 1e+9.
 
Examples:
  plotr '2*T(N/2) + N' 'N*lg(N)'
  plotr -n 1000 -t 10000 'N^2/2' 'T(N-1) + N'
  logplotr 'T(N-1) + T(N-2)' '4; 6; T(N-1) + T(N-2)'
 
USAGE
    exit;
}
 
 
sub longusage {
    print <<USAGE;
 
plotr - plot functions defined by recurrence equations
logplotr - same as plotr, but plots using logarithmic scaling
 
Usage:
  plotr -p -e -n nlimit -t tlimit fn fn ...
  logplotr -p -e -n nlimit -t tlimit fn fn ...
 
Options:
  -e    Print plot as Encapsulated Postscript to plotr.eps if normal exit
  -h    Help.  Prints this explanation of plotr.
  -n    Specifies limit on value of N (x-axis).  Default is 100.
        More than 1000 can be very slow for certain types of recurrences.
  -p    Print plot as Postscript to plotr.ps if normal exit
  -t    Specifies limit on T(N) (y-axis).  Default is 1e+9.
 
Functions:
  Here are a few examples with some explanation of each example.
 
 
  The following plots the recurrence T(N) = 2*T(N/2) + N and the
  function T(N) = N*lg(N)
 
      plotr '2*T(N/2) + N' 'N*lg(N)'
 
  Capital N is used to refer to the input to the function.  In general,
  use single quotes around functions.  Capital T is used to refer to
  the function.  By default, T(0) = 0 and T(1) = 1 (see below to change
  this).  If N/2 is odd, then an average of T((N-1)/2) and T((N+1)/2)
  is used to calculate T(N/2).  To calculate T(N/3), T(N/4), and so on,
  a weighted average is used.
 
 
  The following plots the recurrence T(N) = T(N-1) + N and (N^2)/2 up
  to N = 1000 and T(N) = 10000.
 
      plotr -n 1000 -t 10000 'N^2/2' 'T(N-1) + N'
 
 
  The following plots T(N) = T(N-1) + T(N-2) twice.  Once with the
  default values for T(1) (and T(0)), and a second time using the
  values T(1) = 4 and T(2) = 6.
 
      logplotr 'T(N-1) + T(N-2)' '4; 6; T(N-1) + T(N-2)'
 
  The initial values for the recurrence are separated by semicolons.
 
 
  Note: plotr is implemented in Perl and calls gnuplot for plotting.
  It creates a file for each function that is plotted: tmp0.plt,
  tmp1.plt, and so on.  These are not deleted when plotr exits.
 
  Any Perl arithmetic operation or function can be used, including,
  but not limited to:
 
      cos(x) - cosine of x
      exp(x) - e to the power x, e = 2.718...
      int(x) - integer portion of x, equivalent to floor(x) for x >= 0,
               to get ceiling(x) for x >= 0, use -int(-x)
      log(x) - natural logarithm of x
      sin(x) - sine of x
      sqrt(x) - square root of x
 
  In addition, the following functions are provided.
 
      fact(x) - x!
      lg(x) - logarithm of x, base 2
      ln(x) - natural logarithm of x
      log3(x), log4(x), ... - logarithm of x, base 3, base 4, ...
 
Limitations:
 
  No error checking whatsoever is performed.

  plotr does not solve recurrences.  plotr calculates values recursively,
  and records all intermediate results.  This means that plotr will be
  time- and memory-consuming for large values of n when recurrences
  include terms like T(n-1).  For recurrences that includes terms like
  T(n/b), plotr will generally be quite reasonable for n up to 1,000,000.

Author:
 
  Tom Bylander, 1997.
 
USAGE
    exit;
}
 
 
use Getopt::Std;
if ( getopts('hn:N:t:T:ep') != 1 ) {
    if ( $opt_h ) { longusage; }
    shortusage;
}
 
if ( $opt_h ) { longusage; }
 
$maxfact = 0;
$F[0] = 1;
 
sub fact {
    local($k) = int($_[0]);
    if ( $k > $maxfact ) {
        $F[$k] = &fact($k-1)*$k;
        $maxfact = $k;
    }
    $F[$k];
}
 
sub logb {
    local($base) = $_[0];
    local($x) = $_[1];
    log($x)/log($base);
}

sub recur {
    local($m) = $_[0];
    if ( $m <= 0 ) {
	return 0;
    }
    if ( $m < 1e9 ) {
	local($intm) = int($m);
	if ( $intm != $m ) {
	    return ($m-$intm)*&recur($intm+1) + (1+$intm-$m)*&recur($intm);
	}
    }
    if ( $Tdefined{$m} ) {
	return $T{$m};
    }
    if ( $Trecur{$m} ) {
	return &recur($m-1);
    }
    local($n) = $m;
    $Trecur{$m} = 1;
    $T{$m} = eval($expr);
    delete $Trecur{$m};
    $Tdefined{$m} = 1;
    return $T{$m};
}

sub plotmax {
    local($max) = $_[0];
    local($plotmax) = 10;
    while ( $plotmax < $max ) { $plotmax *= 10; }
    if ( ! $logscale ) {
	if ( $plotmax * 15 / 100 > $max ) { $plotmax = $plotmax * 15 / 100; }
	elsif ( $plotmax / 5 > $max ) { $plotmax /= 5; }
	elsif ( $plotmax * 3 / 10 > $max ) { $plotmax = $plotmax * 3 / 10; }
	elsif ( $plotmax * 2 / 5 > $max ) { $plotmax = $plotmax * 2 / 5; }
	elsif ( $plotmax / 2 > $max ) { $plotmax /= 2; }
	elsif ( $plotmax * 3 / 5 > $max ) { $plotmax = $plotmax * 3 / 5; }
	elsif ( $plotmax * 4 / 5 > $max ) { $plotmax = $plotmax * 4 / 5; }
    }
    elsif ( $plotmax > 1e10 ) {
	if ( $max <= 1e15 ) { $plotmax = 1e15; }
	elsif ( $max <= 1e20 ) { $plotmax = 1e20; }
	elsif ( $max <= 1e30 ) { $plotmax = 1e30; }
	elsif ( $max <= 1e40 ) { $plotmax = 1e40; }
	elsif ( $max <= 1e50 ) { $plotmax = 1e50; }
	elsif ( $max <= 1e60 ) { $plotmax = 1e60; }
	elsif ( $max <= 1e80 ) { $plotmax = 1e80; }
	elsif ( $max <= 1e100 ) { $plotmax = 1e100; }
	elsif ( $max <= 1e150 ) { $plotmax = 1e150; }
	elsif ( $max <= 1e200 ) { $plotmax = 1e200; }
	else { $plotmax = 1e300; }
    }
    if ( $plotmax < 10 ) { $plotmax = 10; }
    return $plotmax;
}
    
 
if ( $opt_n ) {
    $nlimit = $opt_n;
}
elsif ( $opt_N ) {
    $nlimit = $opt_N;
}
else {
    $nlimit = 100;
}
 
if ( $opt_t ) {
    $tlimit = $opt_t;
}
elsif ( $opt_T ) {
    $tlimit = $opt_T;
}
else {
    $tlimit = "1E9";
}

$nmax = 0;
$tmax = 0;
if ( $#ARGV < 0 ) { shortusage; }
foreach $i (0 .. $#ARGV ) {
    $expr = $ARGV[$i];
    if ( $expr =~ /T/ ) {
        $recurrence = 1;
    }
    else {
        $recurrence = 0;
    }
    $expr =~ s/T/\&recur/g;
    $expr =~ s/N/\$n/g;
    $expr =~ s/\^/**/g;
    $expr =~ s/fact/&fact/g;
    $expr =~ s/ln/log/g;
    $expr =~ s/lg/log2/g;
    $expr =~ s/log(\d+)\s*\(/&logb($1,/g;

    @sequence = split(/;/, $expr);
    $expr = pop(@sequence);

    %T = ();
    %Tdefined = ();
    $T{0} = 0;
    $Tdefined{0} = 1;
    if ( $recurrence ) {
	$T{1} = 1;
	$Tdefined{1} = 1;
    }
    $file = "tmp" . $i . ".plt";

    open(TMP, ">$file");
    if ( $#sequence >= 0 ) {
	$n = 1;
	while ( $#sequence >= 0 ) {
	    $T{$n} = shift(@sequence);
	    $Tdefined{$n} = 1;
	    $n++;
        }
    }

    foreach $i (0 .. 100) {
	if ( $i == 0 ) {
	    $n = 1;
	}
	elsif ( $i == 100 ) {
	    $n = $nlimit;
	}
	elsif ( $0 =~ /logplotr/ ) {
	    $n = exp(log($nlimit)*$i/100);
	}
	else {
	    $n = ($nlimit/100)*$i;
	}

	if ( $n <= 1 ) { $n = 1; }
	elsif ( $n < 1e9 ) { $n = int($n); }

	if ( $Tdefined{$n} ) {
	    $val = $T{$n};
	}
	else {
	    $val = eval($expr);
	}

        $T{$n} = $val;
	$Tdefined{$n} = 1;
        printf TMP "%g %g\n", $n, $val;
        if ( $val > $tlimit ) {
            last;
        }
    }
    if ( $val >= $tmax ) {
        $tmax = $val;
    }
    if ( $n >= $nmax ) {
	$nmax = $n;
    }
    close(TMP);
 
}
 
open(GNUPLOT, "| gnuplot");
select((select(GNUPLOT), $| = 1)[0]);
print GNUPLOT "set xlabel \"N\"\n";
print GNUPLOT "set ylabel \"T(N)\"\n";
if ( $0 =~ /logplotr/ ) {
    print GNUPLOT "set logscale\n";
    $logscale = 1;
}

$tmax = &plotmax($tmax);
if ( $tmax >= $tlimit ) {
    $tmax = $tlimit;
}

$nmax = &plotmax($nmax);
if ( $nmax >= $nlimit ) {
    $nmax = $nlimit;
}

if ( $logscale ) {
    printf GNUPLOT "plot [1:%g] [1:%g] ", $nmax, $tmax;
}
else {
    printf GNUPLOT "plot [0:%g] [0:%g] ", $nmax, $tmax;
}

foreach $j (0 .. $#ARGV ) {
    $expr = $ARGV[$j];
    $file = "tmp" . $j . ".plt";
    if ( $j > 0 ) { print GNUPLOT "," }
    print GNUPLOT "\\\n  \"$file\" title \"$expr\" with lines";
}
print GNUPLOT "\n";
 
print "Type return to exit\n";
if ( <STDIN> ) {
    if ( $opt_p ) {
        print GNUPLOT "set term post \"Times-Roman\" 14\n";
        print GNUPLOT "set output \"plotr.ps\"\n";
        print GNUPLOT "replot\n";
        sleep 1;
    }
    if ( $opt_e ) {
        print GNUPLOT "set term post eps \"Times-Roman\" 20\n";
        print GNUPLOT "set output \"plotr.eps\"\n";
        print GNUPLOT "replot\n";
        sleep 1;
    }
    print GNUPLOT "quit\n";
    close(GNUPLOT);
    exit;
}