#!/usr/bin/perl -w

use strict;
use Symbol;
use Scalar::Util;
use clifford;
use Math::Trig qw(pi cot csc asinh acosh atanh);
#
# Must explicitly mention some POSIX functions that are not
# imported by :DEFAULT.
# We need them because otherwise perl would consider them
# built-in, and would not put them in the symbol table,
# so we would not be able to call them by name.
#
# :DEFAULT brings in isatty and some trig functions.
use POSIX qw(:DEFAULT sin cos atan2
              log10 exp);
sub usage {

print <<\EoF;
Desk calculator for Clifford algebra in arbitrarily many
Euclidean dimensions.  (No Minkowski space yet; sorry.)

Usage:  ./cliffer [options]

Command-line options include
  -h      print this message (and exit immediately).
  -v      increase verbosity.
  -i fn   take input from file 'fn'.
  -pre fn take preliminary input from file 'fn'.
  --      take input from STDIN

If no input files are specified with -i or --, the default is an
implicit '--'.  Note that -i and -pre can be used multiple times.
All -pre files are processed before any -i files.

      Advanced usage: If you want to make an input file into a
      self-executing script, you can use "#! /path/to/cliffer -i" as
      the first line.  Similarly, if you want to do some
      initialization and then read from standard input, you can use
      "#! /path/to/cliffer -pre" as the first line.

Ordinary usage example:
  # compound rotation:  two 90 degree rotations
  # makes a 120 degree rotation about the 1,1,1 diagonal:
  echo -e "1 0 0 90° vrml 0 0 1 90° vrml mul @v" | cliffer
Result:
          0.57735 0.57735 0.57735 2.09440 = 120.0000°

Explanation:
  *) Push a rotation operator onto the stack, by
  giving four numbers in VRML format
     X Y Z theta
  followed by the "vrml" keyword.
  *) Push another rotation operator onto the stack,
  in the same way.
  *) Multiply them together using the "mul" keyword.
  *) Pop the result and print it in VRML format using
    the "@v" keyword

On input, we expect all angles to be in radians.  You can convert
  from degrees to radians using the "deg" operator, which can be
  abbreviated to "°" (the degree symbol).  Hint: Alt-0 on some
  keyboards.
  As a special case, on input, a number with suffix "d" (with no
  spaces between the number and the "d") is converted from degrees to
  radians.
    echo "90° sin @" | cliffer
    echo "90 ° sin @" | cliffer
    echo "90d sin @" | cliffer
  are each equivalent to
    echo "pi 2 div sin @" | cliffer

Input words can be spread across as many lines (or as few) as you
  wish.  If input is from an interactive terminal, any error causes
  the rest of the current line to be thrown away, but the program does
  not exit.  In the non-interactive case, any error causes the program
  to exit.

On input, a comma or tab is equivalent to a space.  Multiple spaces
are equivalent to a single space.

Note on VRML format:  X Y Z theta
  [X Y Z] is a vector specifying the axis of rotation,
  and theta specifies the amount of rotation around that axis.
  VRML requires [X Y Z] to be normalized as a unit vector,
  but we are more tolerant;  we will normalize it for you.
  VRML requires theta to be measured in radians.

Also note that on input, the VRML operator accepts either four
  numbers, or one 3-component vector plus one scalar, as in the
  following example.

Same as previous example, with more output:
  echo -e "[1 0 0] 90° vrml dup @v dup @m
           [0 0 1] -90° vrml rev mul dup @v @m" | cliffer
Result:
        1.00000 0.00000 0.00000 1.57080 =  90.0000°
           [  1.00000  0.00000  0.00000 ]
           [  0.00000  0.00000 -1.00000 ]
           [  0.00000  1.00000  0.00000 ]
        0.57735 0.57735 0.57735 2.09440 = 120.0000°
           [  0.00000  0.00000  1.00000 ]
           [  1.00000  0.00000  0.00000 ]
           [  0.00000  1.00000  0.00000 ]

Even fancier:  Multiply two vectors to create a bivector,
then use that to crank a vector:
  echo -e "[ 1 0 0 ] [ 1 1 0 ] mul normalize [ 0 1 0 ] crank @" \
   | ./cliffer
Result:
    [-1, 0, 0]

Another example:  Calculate the angle between two vectors:
  echo -e "[ -1 0 0 ] [ 1 1 0 ] mul normalize rangle @a" | ./cliffer
Result:
    2.35619 = 135.0000°

Example: Powers: Exponentiate a quaternion.  Find rotor that rotates
only half as much:
   echo -e "[ 1 0 0 ] [ 0 1 0 ] mul 2 mul dup rangle @a " \
        " .5 pow dup rangle @a @" | ./cliffer
Result:
    1.57080 =  90.0000°
    0.78540 =  45.0000°
    1 + [0, 0, 1]§

Example:  Take the 4th root using pow, then take the fourth
power using direct multiplication of quaternions:
  echo "[ 1 0 0 ] [ 0 1 0 ] mul dup @v
  .25 pow dup @v dup mul dup mul @v" | ./cliffer
Result
    0.00000 0.00000 1.00000 3.14159 = 180.0000°
    0.00000 0.00000 1.00000 0.78540 =  45.0000°
    0.00000 0.00000 1.00000 3.14159 = 180.0000°

More systematic testing:
   ./cliffer.test1

The following operators have been implemented:
EoF
  listops();
}

######################
## Global variables ##
my $interactive = 0;
my $verbose = 0;
my @stack = ();
my $verb;                       # global, because used in error messages
my @ops = ();
my @private = ();       # same as above, but doesn't show up in help msg
my @lib_unary = ('sin', 'cos', 'tan', 'sec', 'csc', 'cot',
                'sinh', 'cosh', 'tanh',
                'asin', 'acos', 'atan',
                'asinh', 'acosh', 'atanh',
                'ln', 'log2', 'log10', 'exp');
my @lib_binary = ('atan2');
my $deg = pi()/180.;            # one degree, measured in radians
my $asc_hodge = "\247";         # simple ascii § character
my $u2_hodge = "\302$asc_hodge";
my $asc_deg  = "\260";          # simple ascii ° character
my $u2_deg   = "\302$asc_deg";

sub print_vrml {
  my ($A) = @_;
  my ($w, $Ayz, $Azx, $Axy) = @$A;
## convert to the unenlightened normalization that VRML prefers:
  my $norm = sqrt($Ayz*$Ayz + $Azx*$Azx + $Axy*$Axy);
  if ($norm == 0) {
    $Ayz = 1; $norm = 1;
  }
  my $theta = 2.0 * acos($w);
  printf("    %7.5f %7.5f %7.5f %7.5f = %8.4f$u2_deg\n",
     $Ayz/$norm,
     $Azx/$norm,
     $Axy/$norm, $theta, $theta/$deg);
}

## Convert a rotor into a 3x3 rotation matrix

# Print a generic 3x3 matrix
sub print_matrix{
  my ($M) = @_;
  my ($A, $B, $C) = @$M;
  my @rotmat = (@$A, @$B, @$C);

## I hate it when perl formats numbers as "-0.00000".
## We are only printing 5 decimal places, so whacking things
## at the 1e-8 level should be quite harmless:
  for (my $ii = 0;  $ii < $#rotmat; $ii++){
    if (abs($rotmat[$ii]) < 1e-8) {
      $rotmat[$ii] = 0;
    }
  }
  my $f = "      [ %8.5f %8.5f %8.5f ]\n";
  printf("$f$f$f", @rotmat);
}

sub check_depth{
  my ($depth) = @_;
  my $nst = @stack;
  if ($depth > $nst) {
    die "$verb needs $depth things; got only $nst.\n";
  }
}

sub check_number{
  my ($depth) = @_;
  check_depth($depth);
  for (my $ii = 0; $ii < $depth; $ii++){
    my $type = ref($stack[$ii]);
    if ($type ne '') {
      die "$verb needs $depth numbers; got '$type' at depth $ii.\n";
    }
  }
  return 1;
}

# See if the stack contains the required number of things
# of the required type.
sub check_thing{
  my ($depth, $goodtype) = @_;
  check_depth($depth);
  for (my $ii = 0; $ii < $depth; $ii++){
    my $type = myref($stack[$ii]);
    if ($type ne $goodtype) {
      die "$verb needs $depth things of type '$goodtype'; got '$type' at depth $ii.\n";
    }
  }
  return 1;
}

# Call with a clif and the corresponding grademap
sub fmt_vector_row{
    my ($A, $map) = @_;
    my $row = 1;
    my $list = @$map[$row];
    my $vec = [];
    foreach my $slot (@$list){
      push @$vec, @$A[$slot];
    }
    if (codot($vec, $vec) == 0) {
      return '';
    }
    return sprintf('%s%s%s', '[', join(', ', @$vec), ']');
}

# Especially concise printout of clif
# that is already known to have dimensionality D<=3.
sub print_clif3{
  my ($A) = @_;

# all zero is a special case, because of the $didsome elision:
  if (codot($A, $A) == 0) {
    print "    0\n";
    return;
  }
  my $didsome = 0;

  my $compo = @$A;
  my $dim = ilog($compo);
  my @map = grademap($dim);

  my $As = @$A[0];            # the scalar part
  doit: {
    if ($As ne 0) {
      print "    ";
      print $As;
      $didsome++;
    }
    if ($dim == 0){
      last doit;
    }

    my $vec = fmt_vector_row($A, \@map);
    if ($vec ne '') {
      print ($didsome? " + " : "    ");
      print $vec;
      $didsome++;
    }

    if ($dim <= 1) {
      last doit;                # nothing else to do in D=1
    }

# In d=3, print the pseuovector part as the dual of a vector:
    if ($dim == 3) {
      my $vec = fmt_vector_row(hodge($A), \@map);
      if ($vec ne '') {
        print ($didsome? " + " : "    ");
        print $vec, $u2_hodge;
        $didsome++;
      }
    }

# Last but not least, the pseuo-scalar part:
    my $Aps = @$A[$compo-1];
    if ($Aps ne 0) {
      print ($didsome? " + " : "    ");
      print "$Aps$u2_hodge";
      $didsome++;
    }
  }
  print "\n";
}

sub print_clif_rows{
  my ($A) = @_;
  my $compo = @$A;
  my $dim = ilog($compo);
  my @map = grademap($dim);
  for (my $row = 0; $row <= $dim; $row++){
    my $list = $map[$row];
    my $didsome = 0;
    print "    ";
    foreach my $slot (@$list){
      my $val = @$A[$slot];
      my $sign = '+';
      if ($didsome) {
        if ($val < 0) {
          $val *= -1;
          $sign = '-';
        }
        print STDOUT " $sign ";
      }
      print $val, bit_basis($slot);
      $didsome++;
    }
    print "\n";
  }
}

sub print_clif{
  my ($A) = @_;
  if (@$A <= 8) {
    print_clif3($A);
    return;
  }
  print_clif_rows($A);
}

# Format and print an arbitrary item.
# Second argument is a flag ==> also print the ref-type of item.
sub print_general{
  my ($thing, $showtype) = @_;
  my $type = myref($thing);
  if ($showtype) {
    my $dim = -1;
    if ($type eq 'CLIF') {
      $dim = ilog(0+@$thing);
    }
    elsif ($type eq 'VECTOR') {
      $dim = 0+@$thing;
    }
    elsif ($type eq 'QUAT') {
      $dim = 3;
    }
    if ($dim >= 0) {
      printf ("%-7s %-4s:", $type, "D=$dim");
    } else {
      printf ("%-12s:", $type);
    }
  }
  if ($type eq 'ARRAY') {
    my @stuff = @$thing;
    print ("    (", join(', ', @stuff), ")\n");
  }
  elsif ($type eq 'VECTOR') {
    my @stuff = @$thing;
    print ("    [", join(', ', @stuff), "]\n");
  }
  elsif ($type eq '--MARK--'){
    print "    [\n";
  }
  else {
    $thing = promote($thing);
    if (@$thing > 8) {
      print "\n";
    }
    print_clif($thing);
  }
}

######################################################################
## populate the list of operators ............

sub listops{
  foreach my $cmd (@ops) {
    my $type = ref($cmd);
    if (!$type){
      printf "=== $cmd\n";
    }
    else {
      printf "%-7s %s\n", @$cmd[0], @$cmd[1];
    }
  }
  print "=== Math library functions:\n";
  my $didsome = 0;
  foreach my $cmd (@lib_unary, @lib_binary){
    if (($didsome % 6) == 0) {
      if ($didsome) {
        print "\n";
      }
      $didsome = 0;
      printf "%7s", "";
    }
    printf "%-6s ", $cmd;
    $didsome++;
  }
  if ($didsome) {
    print "\n";
  }
}

push @private, ['?', "help message", \&usage];
push @private, ['usage', "help message", \&usage];
push @ops, ['help', "help message", \&usage];
push @ops, ['listops', "list all operators", \&listops];

push @ops, "Unary operators";

push @ops, ['pop', "remove top item from stack", sub {
  check_depth(1);
  shift @stack;
}];

push @ops, ['neg', "negate: multiply by -1", sub {
  check_number(1);
  $stack[0] *= -1;
}];

push @ops, ['deg', "convert number from radians to degrees", sub {
  check_number(1);
  $stack[0] *= &pi/180;
}];

push @ops, ['dup', "duplicate top item on stack", sub {
  check_depth(1);
  unshift @stack, $stack[0];
}];

push @ops, ['gorm', "gorm i.e. scalar part of V~ V", sub {
  check_depth(1);
  $stack[0] = codot($stack[0], $stack[0]);
}];

push @ops, ['norm', "norm i.e. sqrt(gorm}", sub {
  check_depth(1);
  $stack[0] = norm($stack[0]);
}];

push @ops, ['normalize', "divide top item by its norm", sub {
  check_depth(1);
  $stack[0] = normalize($stack[0]);
}];

push @ops, ['rev', "clifford '~' operator, reverse basis vectors", sub {
  check_depth(1);
  $stack[0] = rev($stack[0]);
}];

sub do_hodge {
  check_depth(1);
  $stack[0] = hodge($stack[0]);
}

push @ops, ['hodge',
  "hodge dual aka unary '$u2_hodge' operator; alt-' on some keyboards",
sub {do_hodge}
];

push @ops, ['gradesel', "given C and s, find the grade-s part of C", sub{
  check_depth(2);
  check_number(1);
  my $B = shift @stack;
  check_thing(1, 'CLIF');
  my $A = $stack[0];
  $stack[0] = gradesel($A, $B);
}];

push @ops, ['rangle', "calculate rotor angle", sub {
  check_thing(1, 'QUAT');
  $stack[0] = rangle($stack[0]);
}];

push @ops, "Binary operators";

push @ops, ['exch', "exchange top two items on stack", sub {
  check_depth(2);
  ($stack[0], $stack[1]) = ($stack[1], $stack[0]);
}];

push @ops, ['codot', "multiply corresponding components, then sum", sub {
  check_depth(2);
  my $B = shift @stack;
  $stack[0] = codot($stack[0], $B);
}];

push @ops, ['add', "add top two items on stack", sub {
  check_depth(2);
  my $B = shift @stack;
  $stack[0] = add($stack[0], $B);
}];

push @ops, ['sub', "sub top two items on stack", sub {
  check_depth(2);
  my $B = shift @stack;
  $stack[0] = add($stack[0], mul_by_scalar($B, -1));
}];

push @ops, ['mul',
"multiply top two items on stack (in subspace if possible)", sub {
  check_depth(2);
  my $B = shift @stack;
  my $A = $stack[0];
  my $Atype = myref($A);
  my $Btype = myref($B);
  if (0) {}
  elsif ($Atype eq 'plain_number') {
    $stack[0] = mul_by_scalar($B, $A);
  }
  elsif ($Btype eq 'plain_number') {
    $stack[0] = mul_by_scalar($A, $B);
  }
  elsif ($Atype eq 'VECTOR' && $Btype eq 'VECTOR'
        && (@$A == 3) && (@$B == 3) ) {
    $stack[0] = vec3mul($A, $B);
  }
  elsif ($Atype eq 'QUAT' && $Btype eq 'QUAT') {
    $stack[0] = quatmul($A, $B);
  }
  else {
# If all else has failed, promote them both and
# invoke the general case:
    $stack[0] = clifmul(promote($A), promote($B), 0);
  }
}];

push @ops, ['cmul', "promote A and B to clifs, then multiply them", sub{
  check_depth(2);
  my $B = shift @stack;
  my $A = $stack[0];
  $stack[0] = clifmul(promote($A), promote($B), 0);
}];

push @ops, ['div', "divide clif A by scalar B", sub{
  check_depth(2);
  check_number(1);
  my $B = shift @stack;
  my $A = $stack[0];
    $stack[0] = mul_by_scalar($A, 1/$B);
}];

push @ops, ['dot', "promote A and B to clifs, then take dot product", sub{
  check_depth(2);
  my $B = shift @stack;
  my $A = $stack[0];
  $stack[0] = clifmul(promote($A), promote($B), -1);
}];

sub do_wedge {
  check_depth(2);
  my $B = shift @stack;
  my $A = $stack[0];
  $stack[0] = clifmul(promote($A), promote($B), 1);
}

push @ops, ['wedge', "promote A and B to clifs, then take wedge product",
sub{do_wedge}];

push @ops, ['cross', "the hodge of the wedge (familiar as cross product in 3D)",
  sub{do_wedge; do_hodge}
];

push @ops, ['crank', "calculate R~ V R", sub {
  check_thing(1, 'VECTOR');
  my $V = shift @stack;
  check_thing(1, 'QUAT');
  $stack[0] = crank($stack[0], $V);
}];

sub do_pow{
  check_number(1);
  my $S = shift @stack;
  my $X = $stack[0];
  my $Stype = myref($S);
  my $Xtype = myref($X);
  if ($Xtype eq 'QUAT') {
    $stack[0] = quat_pow($X, $S);
  }
  elsif ($Xtype eq 'plain_number') {
    $stack[0] = $X ** $S;
  }
  else {
    die "Sorry, don't know how to take powers using '$Xtype' and '$Stype'\n";
  }
}

push @ops, ['pow', "calculate Nth power of scalar or quat",
  sub{do_pow}
];

push @ops, ['sqrt', "calculate square root of power of scalar or quat",
  sub{
    unshift @stack, 0.5;
    do_pow;
  }
];

push @ops, "Constructors";

push @ops, ['[', "mark the beginning of a vector", sub{
    unshift @stack, [];
    bless $stack[0], '--MARK--';
}];

push @ops, [']', "construct vector by popping to mark", sub{
  my $rslt = [];
  popper: for (;;){
    if (!@stack) {
      die "constructor ']' could not find mark '['\n";
    }
    if (ref($stack[0]) eq '--MARK--') {
      shift @stack;
      last popper;
    }
    unshift @$rslt, shift @stack;
  }
  bless $rslt, 'VECTOR';
  unshift @stack, $rslt;
}];

push @ops, ['unpack',
"unpack a vector, quat, or clif; push its contents (normal order)",
sub{do_unpack()}
];

sub do_unpack{
  check_depth(1);
  my $A = shift @stack;
  my $type = myref($A);
  if ($type eq 'plain_number'){
    unshift @stack, $A;         # just put it back
  }
  elsif ($type eq 'QUAT' || $type eq 'VECTOR') {
    foreach my $comp (@$A) {
      unshift @stack, $comp;
    }
  }
  elsif ($type eq 'CLIF'){
# code cloned from print_clif_rows:
    my $compo = @$A;
    my $dim = ilog($compo);
    my @map = grademap($dim);
    for (my $row = 0; $row <= $dim; $row++){
      my $list = $map[$row];
      foreach my $slot (@$list){
        my $val = @$A[$slot];
        unshift @stack, $val;
      }
    }
  }
  else {
    die "Don't know how to unpack item of type '$type'\n";
  }
}

push @private, ['bunpack',
"unpack a vector, quat, or clif; push its contents (bitcode order)",
sub{
  check_depth(1);
  my $A = shift @stack;
  my $type = myref($A);
  if ($type eq 'plain_number'){
    unshift @stack, $A;         # just put it back
  }
  elsif ($type eq 'QUAT'
    || $type eq 'VECTOR'
    || $type eq 'CLIF'
    ) {
    foreach my $comp (@$A) {
      unshift @stack, $comp;
    }
  }
  else {
    die "Don't know how to bunpack item of type '$type'\n";
  }
}];

push @ops, ['dimset', "project object onto N-dimensional Clifford algebra",
sub{
  check_number(1);
  my $dim = shift @stack;
  check_depth(1);
  my $A = promote(shift @stack);
  my $ncomp = 2**$dim;        # number of components in desired result
  my $todo = min($ncomp, 0+@$A);
  my $rslt = [];
  for (my $ii = 0; $ii < $todo; $ii++){
    @$rslt[$ii] = @$A[$ii];
  }
  for (my $ii = $todo; $ii < $ncomp; $ii++){
    @$rslt[$ii] = 0;
  }
  bless $rslt, 'CLIF';
  unshift @stack, $rslt;
}];

push @ops, ['unbave', "top unit basis vector in N dimensions", sub{
  check_number(1);
  my $dimen = $stack[0];
  if ($dimen <= 0) {
    die "No vectors in dimension $dimen.\n";
  }
  my $ndx = $dimen-1;   # x0 is a vector in 1 dimension
                        # x0 and x1 are vectors in 2 dimensions
  my $rslt = [];
  my $compo = 2**$dimen;
  for (my $ii = 0; $ii < $compo; $ii++){
    @$rslt[$ii] = 0;
  }
  @$rslt[2**$ndx] = 1;
  bless $rslt, 'CLIF';
  $stack[0] = $rslt;
}];

push @ops, ['ups', "unit pseudo-scalar in N dimensions", sub{
  check_number(1);
  my $dimen = $stack[0];
  if ($dimen <=0) {
    die "No vectors in dimension $dimen.\n";
  }
  my $rslt = [];
  my $compo = 2**$dimen;
  for (my $ii = 0; $ii < $compo; $ii++){
    @$rslt[$ii] = 0;
  }
  @$rslt[$compo-1] = 1;
  bless $rslt, 'CLIF';
  $stack[0] = $rslt;
}];

push @ops, ['pi', "push pi onto the stack", sub{
  unshift @stack, pi();
}];

push @ops, ['vrml', "construct a quaternion from VRML representation x,y,z,theta", sub {

## Let the VRML representation be X, Y, Z, theta,
## where [X,Y,Z] is a unit vector.
## Here (Aw; Ax, Ay, Az) is the multivector representation,
##  i.e. quaternion representation, *not* the VRML representation.
##  Specifically, w1=cos(theta/2) and x1 is X*sin(theta/2).
##  Note that the whole multivector (w1; x1, y1, z1)
##  (including the scalar part w1, not just the bivector part)
##  is normalized to unity.

  check_number(1);
  my $theta = shift @stack;
  check_depth(1);
  my $type = myref($stack[0]);
  if ($type eq 'VECTOR' && 0+@{$stack[0]} == 3) {
    do_unpack();
  }
  check_number(3);
  my $Az = shift @stack;
  my $Ay = shift @stack;
  my $Ax = shift @stack;
# should already be normalized, but we force the issue:
  my $norm = sqrt($Ax*$Ax + $Ay*$Ay + $Az*$Az);
  if ($norm == 0) {
    $Ax = 1; $norm = 1; $theta = 0;
  }
  my $ns = sin($theta/2.0) / $norm;       # normalized sine
  my $rotor = [cos($theta/2.0), $Ax*$ns, $Ay*$ns, $Az*$ns];
  bless $rotor, 'QUAT';
  unshift @stack, $rotor;
}];

push @ops, ['clif', "take a vector in D=2**n, construct a clif in D=n\n"
. "\tNote:  You can do the opposite via '[ exch unpack ]'",
sub {
    check_thing(1, 'VECTOR');
    my $A = shift @stack;
    my $rslt = [];
# code cloned from print_clif_rows:
    my $compo = @$A;
    my $dim = ilog($compo);
    my @map = grademap($dim);
    for (my $row = 0; $row <= $dim; $row++){
      my $list = $map[$row];
      foreach my $slot (@$list){
        my $val = @$A[$slot];
        push @$rslt, $val;
      }
    }
    bless $rslt, 'CLIF';
    unshift @stack, $rslt;
}];

push @ops, "Printout operators";

my $basis_mode = 0;

push @ops, ['setbasis', "set basis mode, 0=abcdef 1=xyzabc", sub{
  check_number(1);
  $basis_mode = shift @stack;
}];

sub bit_basis{
  my ($code, $basis) = @_;
  my $rslt = '';
  if (! defined $basis) {
    $basis = 'abcdefghijklmnopqrstuvwxyz';
    if ($basis_mode) {
      $basis = 'xyz' . $basis;
    }
  }
  my @list = bitcode_unpack($code);
  foreach my $bitpos (@list) {
    $rslt .= substr($basis, $bitpos, 1);
  }
  return $rslt;
}

push @ops, ['dump', "show everything on stack, leave it unchanged", sub {
  my $nst = @stack;
  ## print "$nst items:\n";
  for (my $ii = 0; $ii < $nst; $ii++) {
    printf "%2d: ", $ii;
    print_general($stack[$ii], 1);
  }
}];

push @ops, ['@', "compactly show item of any type, D=3 (then remove it)", sub {
  check_depth(1);
  print_general (shift @stack);
}];

push @ops, ['@m', "show quaternion, formatted as a rotation matrix (then remove it)", sub {
  check_thing(1, 'QUAT');
  print_matrix(rotmat(shift @stack));
}];

push @ops, ['@v', "show quaternion, formatted in VRML style (then remove it)", sub {
  check_thing(1, 'QUAT');
  print_vrml (shift @stack);
}];

push @ops, ['@a', "show angle, formatted in radian and degrees (then remove it)", sub {
  check_number(1);
  my $theta = shift @stack;
  printf("    %7.5f = %8.4f$u2_deg\n", $theta, $theta/$deg);
}];

push @ops, ['@x', "print clif of any grade, row by row", sub{
  check_depth(1);
  print_clif_rows (promote(shift @stack));
}];

sub tabulate_grade_of_product{
  my $dim = 3;
  my $comp = 2**$dim;
  print "grade of product:\n";
# Print column headers
  printf "%9s", '';
  for (my $jj = 0; $jj < $comp; $jj++) {
    printf " %3d     ", $jj;
  }
  print "\n";
  printf "%9s", '(';
  for (my $jj = 0; $jj < $comp; $jj++) {
    printf " %3d     ", bladegrade($jj);
  }
  print ")\n\n";

  for (my $ii = 0; $ii < $comp; $ii++) {
    my $iigrade = bladegrade($ii);
    printf " %3d(%2d):", $ii, $iigrade;
    for (my $jj = 0; $jj < $comp; $jj++) {
      my $jjgrade = bladegrade($jj);
      my $prod = $ii ^ $jj;
      my $flags = '';
      my $pgrade = bladegrade($prod);
      $flags .= ($pgrade == abs($iigrade - $jjgrade)) ? 'd' : '_';
      $flags .= ($pgrade == abs($iigrade + $jjgrade)) ? 'w' : '_';
      $flags .= (sign_of_product($ii, $jj) > 0) ? '+' : '-';
#     $flags .= (sign_of_product($jj, $ii) > 0) ? '+' : '-';
      printf " %3d %-4s", $pgrade, $flags;
    }
    print "\n";
  }
}

sub test_blademap{
  my $dim = 3;          # for testing
  my $comp = 2**$dim;
  print "blademap:\n";
  for (my $ii = 0; $ii < $comp; $ii++) {
    printf "%3d", $ii;
    my @map = bladeMap($ii, $dim);
    foreach my $jj (@map){
      printf " %3d", $jj;
    }
    print "\n";
  }
}

push @private, ['play', "experiment with counters, general clif multiply", sub{
  print "Sign of product: ",
     sign_of_product(1, 2), ", ",
     sign_of_product(2, 1), "\n";
  test_blademap();
  print "...\n";
  tabulate_grade_of_product();
}];

push @private, ['perm', "experiment with permutations", sub{
  my $dim = 4;
  my $universe = 2**$dim - 1;
  my $code1 = 2;
  my $code2 = $universe ^ $code1;
  my @list1 = bitcode_unpack($code1);
  my @list2 = bitcode_unpack($code2);
  print "checking: ",         join(',', @list1, @list2),
        " == ",  sign_of_perm([@list1, @list2]), "\n";
}];


######################################################################

sub log2{
  return log(shift @_) / log(2);
}


sub ln{
  return log(shift @_);
}

###
# Remap a few special symbols.
# See also suffix parse (below).
#
my %special_syms = (
 'VRML'         => 'vrml',
 $u2_hodge      => 'hodge',
 $u2_deg        => 'deg',
 "\317\200"     => 'pi',        # unicode π aka pi => pi
 "\342\213\200" => 'wedge',     # unicode ⋀ aka xwedge => wedge
 "\342\213\205" => 'dot',       # unicode ⋅ aka sdot => dot
 "\302\267"     => 'dot'        # unicode · aka middot => dot
);

main: {
  my @infiles = ();
  my @prefiles = ();

  my @args = @ARGV;
  while (@args) {
    my $arg = shift @args;
    if ($arg eq '-h'){
      usage();
      exit 0;
    } elsif ($arg eq '-v'){
      $verbose++;
    } elsif ($arg eq '--'){
      push @infiles, $arg;
    } elsif ($arg eq '-i'){
      if (!@args) {
        die "Option '-i' requires an argument\n";
      }
      push @infiles, shift @args;
    } elsif ($arg eq '-pre'){
      if (!@args) {
        die "Option '-pre' requires an argument\n";
      }
      push @prefiles, shift @args;
    } else {
      die "Unrecognized option '$arg';  for help try cliffer -h\n";
    }
  }
  if (!@infiles) {
    push @infiles, '--';
  }
  @stack = ();          # start with empty stack
  my $inch = Symbol::gensym;

  filer: foreach my $infile (@prefiles, @infiles) {
   if ($infile eq '--') {
     open($inch, '<&', 'STDIN')
       || die "Could not dup STDIN.\n";
   } else {
     open($inch, '<', $infile)
       || die "Could not open input file '$infile'.\n";
   }
   $interactive = isatty($inch) || 0;

   liner: for (;;) {
    if ($interactive) {
      print STDERR "> ";
    }
    if (eof($inch)) {
      last liner;
    }
    my $line = <$inch>;
    chomp $line;
    my $word;
    my $lookahead = '';
# Suffix parse:
# ?? There are probably cleverer ways to parse commas and
# square brackets using split() with "collector" patterns.
# See ../computer/chem-parse.pl
# See also %special_syms.
    $line =~ s/$u2_deg/ deg/g;
    $line =~ s/$u2_hodge/ hodge/g;
    $line =~ s/$asc_deg/ deg/g;
    $line =~ s/$asc_hodge/ deg/g;
    $line =~ s'[[]'[ 'g;
    $line =~ s']' ]'g;
    $line =~ s',' 'g;
    my @words = split(' ', $line);
    worder: for (;;) {
      my $degree = 0;
      if ($lookahead ne '') {
        $word = $lookahead;
        $lookahead = '';
      } else {
        if (!@words) {
          last worder;
        }
        $word = shift @words;
        if ($word eq '') {
          next worder;
        }
# strip comments:
        if ($word =~ m'^#'){
          next liner;
        }
      }
# Note $verb is a global variable, useful in debug messages:
      $verb = $word;

# Remap a few special symbols:
      if (exists $special_syms{$word}) {
        $verb = $word = $special_syms{$word};
      }

# strip out degree symbol (alt-zero)
# do this BEFORE seeing if $word looks_like_number:
      if ($word =~ s/[d]$// ){
         $degree++;
      }
      if (Scalar::Util::looks_like_number($word)){
        if ($degree) {
          $word *= &pi/180;
        }
        unshift @stack, $word;
        next worder;
      }
# Here if the word happens to be a verb.

      my $found_cmd = '';
      finder: {
        foreach my $cmd (@ops, @private) {
          my $type = ref($cmd);
          if (!$type){
            # skip sectioning comments
          }
          elsif ($verb eq @$cmd[0]){
            $found_cmd = @$cmd[2];
            last finder;
          }
        }
        foreach my $cmd (@lib_unary) {
          if ($verb eq $cmd) {
              $found_cmd = \&unary_wrapper;
              last finder;
          }
        }
        foreach my $cmd (@lib_binary) {
          if ($verb eq $cmd) {
              $found_cmd = \&binary_wrapper;
              last finder;
          }
        }
        print STDERR "Unrecognized verb '$verb'; try 'help' or 'listops'\n";
        if (!$interactive){
          exit 1;
        }
        next liner;
      }

      eval {
        &$found_cmd($verb);
        1;                  # always true unless there is an exception
      } or do {             # catch any exception
        msger: for my $msg (split("\n", $@)) {
# Remove nl (if any) then add nl,
# so this works whether or not the msg had a nl when it was thrown.
          chomp $msg;
# We don't want to see line numbers
# in msgs generated by library functions:
          if ($msg =~ m'^Died at') {
            next msger;
          }
          print STDERR "$msg\n";
        }
        if (!$interactive){
          exit 1;
        }
        next liner;
      };
    }  ## end loop over words
   }   ## end loop over lines in file
   close $inch;
  }     ## end loop over input files
  if ($interactive){            # erase the prompt for absent input
    print STDERR  "\r  \r";
  }
  exit 0;
}

############
# Some simple wrappers, for applying library functions to our stack:

sub check_result_num{
  my ($cmd) = @_;
  my $type = ref($stack[0]);
  if ($type ne ''){
    die "Hmmm, '$cmd' returned a result of type '$type';\n"
      . "should have been a plain number.\n";
  }
}

sub unary_wrapper{
  check_number(1);
  my ($cmd) = @_;
 ## one step of indirection, to get around "strict refs":
  $stack[0] = &{\&$cmd}($stack[0]);
  check_result_num($cmd);
}

sub binary_wrapper{
  check_number(2);
  my ($cmd) = @_;
  my $second = shift @stack;
 ## one step of indirection, to get around "strict refs":
  $stack[0] = &{\&$cmd}($stack[0], $second);
  check_result_num($cmd);
}