package clifford;

use Exporter;
@ISA = ('Exporter');
@EXPORT = ('rotmat', 'quat_pow', 'crank', 'rangle',
        'hodge', 'sign_of_perm',
        'vector_matrix_prod', 'codot',
        'myref', 'rev', 'rev_clif', 'rev_quat',
        'norm', 'normalize',
        'mul_by_scalar', 'promote',
        'vec3mul', 'quatmul', 'clifmul', 'sign_of_product',
        'add', 'copy', 'gradesel',
        'bladeMap', 'grademap', 'bladegrade',
        'bitcode_unpack',
        'min', 'ilog', 'plog_p1');

use Math::Trig;

# pseudo-logarithm plus 1, i.e. minimum number of bits
# rquired to represent a number
#     plog_p1(0) = 0
#     plog_p1(1) = 1
#     plog_p1(4) = plog_p1(7) = 3
#     plog_p1(8) = plog_p1(15) = 4
#
sub plog_p1{
  my ($arg) = @_;
  my $rslt = 0;
  while ($arg) {
    $rslt++;
    $arg >>= 1;
  }
  return $rslt;
}

# requires arg to be a power of 2;
# returns the base-2 logarithm, always an integer
sub ilog{
  my ($arg) = @_;
  my $rslt = plog_p1($arg-1);
  if ($arg != 2**$rslt) {
    die "ilog: arg should be a power of 2, not '$arg'\n";
  }
  return $rslt;
}

sub myref{
  my ($thing) = @_;
  my $type = ref($thing);
  if ($type eq '') {
    $type = 'plain_number';
  }
  return $type;
}

sub rotmat {
  my ($A) = @_;
  my ($w, $Ayz, $Azx, $Axy) = @$A;
  my $theta = 2.0 * acos($w);

  my $rotmat = [
 [$w*$w+$Ayz*$Ayz - $Azx*$Azx - $Axy*$Axy, -2*($w*$Axy - $Ayz*$Azx)           ,  2*($w*$Azx + $Axy*$Ayz)],
 [ 2*($w*$Axy + $Ayz*$Azx)       , $w*$w+$Azx*$Azx - $Axy*$Axy - $Ayz*$Ayz    , -2*($w*$Ayz - $Azx*$Axy)],
 [-2*($w*$Azx - $Axy*$Ayz)       ,     2*($w*$Ayz + $Azx*$Axy), $w*$w+$Axy*$Axy - $Azx*$Azx - $Ayz*$Ayz]];

  bless $rotmat, 'MATRIX';
  return $rotmat;
}

sub quat_pow{
  my ($Q, $pow) = @_;
  my $Qnorm = norm($Q);
  my $rslt = [1, 0, 0, 0];
  bless $rslt, 'QUAT';
  if ($pow == 0) {
    return $rslt;
  }
  my ($As, @Ab) = @{$Q};
  my $bnorm = norm([@Ab]);
  my $theta = atan2($bnorm, $As);
  my $Abunit;
  if ($bnorm > 0) {
    $Abunit = mul_by_scalar([@Ab], 1/$bnorm)
  } else {
    $Abunit = [1, 0, 0];
    # nonfatal warning:
    print STDERR "Plane of rotation undefined; assuming XY plane\n";
  }
  $theta *= $pow;
  #xx print STDERR "theta: $theta  bnorm: $bnorm  abunit: @$Abunit[0] @$Abunit[1]\n";
  my $Bb = mul_by_scalar($Abunit, sin($theta));
  $rslt = [cos($theta), @$Bb];
  bless $rslt, 'QUAT';
  #xx print STDERR "rslt: ", join(",, ", @$rslt), "\n";
  my $foo = mul_by_scalar($rslt, $Qnorm ** $pow);
  #xx print STDERR "foo: ", join(",, ", @$foo), "\n";
  return $foo;
}

sub crank{
  my ($A, $B) = @_;
  return matrix_vector_prod(rotmat($A), $B);
}

# The rotor angle of quaternion.
sub rangle{
  my ($Q) = @_;
  my ($As, @Ab) = @{$Q};
  my $bnorm = norm([@Ab]);
  return atan2($bnorm, $As);
}

sub xxxhodge{
  my ($A) = @_;
  $A = promote($A);
  my $comp = @$A;
  my $dim = ilog($comp);
  if ($dim != 3) {
    die "hodge only implemented for 3 dimensions not $dim.\n";
  }

# New clif packing rule
# 0, 1, 2, 3,  4, 5!,  6,  7
# .  x, y, xy, z, xz, yz, xyz

  my ($As, $Ax, $Ay, $Axy, $Az, $Axz, $Ayz, $Aps) = @$A;

  my $rslt = [$Aps, $Ayz, -$Axz, $Az,
                    $Axy,  $Ay,  $Ax,  $As];
  bless $rslt, 'CLIF';
  return $rslt;
}

sub sign_of_perm{
  my ($arg) = @_;
  my @list = @$arg;
  my $rslt = 1;
  my $dim = @list;
  for (my $ii = 0; $ii < $dim; $ii++){
    my $remaining = @list;
    #xx print STDERR "sign_of_perm... dim: $dim  remaining: $remaining  ",
    #xx     join(',', @list), "\n";
    for (my $jj =0; $jj < $remaining; $jj++){
      if ($list[$jj] == $ii) {
        splice @list, $jj, 1;
        goto done;
      }
      $rslt *= -1;
    }
    die "Missing item '$ii' from permutation in sign_of_perm.\n";
done:
  }
  return $rslt;
}

# Returns a list of numbers indicating
# which bits are set in the bitcode.
# Numbers are in normal order.
sub bitcode_unpack{
  my ($code) = @_;
  my @rslt = ();
  for (my $ii = 0; $code; $ii++){
    if ($code & 1) {
      push @rslt, $ii;
    }
    $code >>= 1;
  }
  return @rslt;
}

sub gradesel{
  my ($A, $goodgrade) = @_;
  $A = promote($A);
  my $comp = @$A;
  my $dim = ilog($comp);
  my @map = grademap($dim);
  my $rslt = [];
  for (my $row = 0; $row <= $dim; $row++){
    my $list = $map[$row];
    foreach my $code1 (@$list){
      @$rslt[$code1] = ($row == $goodgrade) ?
        @$A[$code1] : 0;
    }
  }
  bless $rslt, 'CLIF';
  return $rslt;
}

# the outer part of this code runs somewhat parallel to rev_clif
sub hodge{
  my ($A) = @_;
  $A = promote($A);
  my $comp = @$A;
  my $dim = ilog($comp);
  my $univ = 2**$dim - 1;
  my @map = grademap($dim);
  my $rslt = [];
  for (my $row = 0; $row <= $dim; $row++){
    my $list = $map[$row];
    foreach my $code1 (@$list){
      my $code2 = $univ ^ $code1;
      my @list1 = bitcode_unpack($code1);
      my @list2 = bitcode_unpack($code2);
      @$rslt[$code1] = sign_of_perm([@list1, @list2]) * @$A[$code2];
    }
  }
  bless $rslt, 'CLIF';
  return $rslt;
}

sub matrix_vector_prod{
  my ($M, $V) = @_;
  my ($row1, $row2, $row3) = @$M;
  my $rslt = [
    codot($row1, $V),
    codot($row2, $V),
    codot($row3, $V),
  ];
  bless $rslt, 'VECTOR';
  return $rslt;
}

sub min{
  my @arg = @_;
  if (!@arg) {
    die "min() requires at least one argument.\n";
  }
  my $rslt = shift @arg;
  while (@arg){
    my $x = shift @arg;
    if ($x < $rslt) {
      $rslt = $x;
    }
  }
  return $rslt;
}

# multiply corresponding components,
# which is *not* the Clifford dot product
sub codot{
  my ($A, $B) = @_;
  my $tA = myref($A);
  my $tB = myref($B);
  checkit: {
    if ($tA eq $tB){
      last checkit;
    }
    if ($tA eq 'ARRAY' && $tB eq 'VECTOR') {
      last checkit;
    }
    if ($tA eq 'VECTOR' && $tB eq 'ARRAY') {
      last checkit;
    }
    die "Codot: mismatched types: '$tA' versus '$tB'\n";
  }
  my $ncomp = min(0+@$A, 0+@$B);            # number of components
  my $codot = 0;
  for (my $ii = 0; $ii < $ncomp; $ii++){
    $codot += @$A[$ii] * @$B[$ii];
  }
  return $codot;
}

# we require $B to be a scalar,
# but $A can be almost anything,
# either a scalar or any kind of array of components
sub mul_by_scalar{
  my ($A, $B) = @_;
  my $type = myref($A);
  if ($type eq 'plain_number'){
    return $A * $B;
  }
  my $ncomp = @$A;
  my $C = [];
  for (my $ii = 0; $ii < $ncomp; $ii++){
    @$C[$ii] = @$A[$ii] * $B;
  }
  bless $C, ref($A);
  return $C;
}

# calculate the norm of an array-object such as a quaternion
sub norm{
  my ($A) = @_;
  return sqrt(codot($A, $A));
}

# return a normalized version of an array-object such as a quaternion
sub normalize{
  my ($A) = @_;
  return mul_by_scalar($A, 1/norm($A));
}

sub bit_flip{
  my ($code, $dimen) = @_;
  my $rslt = 0;
  for (my $ii = 0; $ii < $dimen; $ii++){
    $rslt <<= 1;
    if ($code & 1) {
      $rslt |= 1;
    }
    $code >>= 1;
  }
  return $rslt;
}

my $Xdimen = 4;

sub flip_comp{
  return bit_flip($b, $Xdimen) <=> bit_flip($a, $Xdimen);
}

sub grademap{
  my ($dimen) = @_;
  $Xdimen = $dimen;
  my $compo = 2**$dimen;
  my @rslt = ();
  for (my $dim = 0; $dim <= $dimen; $dim++){
    $rslt[$dim] = [];
  }
  for (my $ii = 0; $ii < $compo; $ii++){
    my $grade = bladegrade($ii);
    push @{$rslt[$grade]}, $ii;
  }
  for (my $dim = 0; $dim <= $dimen; $dim++){
    @{$rslt[$dim]} = sort {flip_comp} @{$rslt[$dim]};
  }
  return @rslt;
}

# reverse of a clif.
sub rev_clif{
  my ($A) = @_;
  my $comp = @$A;
  my $dim = ilog($comp);
  my @map = grademap($dim);
  my $rslt = [];
  for (my $row = 0; $row <= $dim; $row++){
    my $sign = (int($row/2) & 1) ? -1 : 1;
    my $list = $map[$row];
    my $didsome = 0;
    foreach my $slot (@$list){
      @$rslt[$slot] = $sign * @$A[$slot];
    }
  }
  bless $rslt, 'CLIF';
  return $rslt;
}

# reverse of a quaternion
sub rev_quat{
  my ($A) = @_;
  my ($As, $Ayz, $Azx, $Axy) = @$A;
  my $rslt = [$As, -$Ayz, -$Azx, -$Axy];
  bless $rslt, 'QUAT';
  return $rslt;
}

# this version works for either clifs or quats:
sub rev{
  my ($A) = @_;
  if (ref($A) eq 'QUAT') {
    return rev_quat($A);
  }
  return rev_clif(promote($A));
}

# Multiply two three-dimensional vectors, to form a QUAT:
sub vec3mul {
  my ($A, $B) = @_;
  my ($Ax, $Ay, $Az) = @$A;
  my ($Bx, $By, $Bz) = @$B;
  my $rslt = [
    $Ax*$Bx + $Ay*$By + $Az*$Bz,
    $Ay*$Bz - $Az*$By,
    $Az*$Bx - $Ax*$Bz,
    $Ax*$By - $Ay*$Bx
  ];
  bless $rslt, 'QUAT';
  return $rslt;
}

# Multiply two things of type QUAT.
# This includes (but is not limited to) rotors.
# Result is a QUAT.
sub quatmul {
  my ($A, $B) = @_;
  my ($Aw, $Ayz, $Azx, $Axy) = @$A;
  my ($Bw, $Byz, $Bzx, $Bxy) = @$B;
  my $rslt = [
    $Aw*$Bw - $Ayz*$Byz - $Azx*$Bzx - $Axy*$Bxy,
    $Aw*$Byz + $Bw*$Ayz + ($Axy*$Bzx - $Azx*$Bxy),
    $Aw*$Bzx + $Bw*$Azx + ($Ayz*$Bxy - $Axy*$Byz),
    $Aw*$Bxy + $Bw*$Axy + ($Azx*$Byz - $Ayz*$Bzx)
  ];
  bless $rslt, 'QUAT';
  return $rslt;
}

sub sign_of_product{
  my ($ii, $jj) = @_;
  my $dim = plog_p1($ii | $jj);
  my @map = bladeMap($jj, $dim);
  my $sign = 1;
  for (my $kk = 0; $kk < $dim; $kk++) {
    my $mask = 2**$kk;
    if ($ii & $mask){
      if ($map[$kk] & 1){
        $sign *= -1;
      }
    }
  }
  return $sign;
}

# Call with $blade using bitcode representation
# We just count the bits that are set in it.
sub bladegrade{
  my ($blade) = @_;
  my $runtot = 0;
  while ($blade) {
    if ($blade & 1) {
      $runtot++;
    }
    $blade >>= 1;
  }
  return $runtot;
}

# The map has 1+$dim elements, from 0 through $dim inclusive.
# Each entry shows the number of factors to the "right" of here.
# Therefore $map[$dim] is the grade of the blade, i.e. the
# total number of factors.
sub bladeMap{
  my ($blade, $dim) = @_;
  my @map = ();
  my $runtot = 0;
  for (my $jj = 0; $jj < $dim; $jj++){
    $map[$jj] = $runtot;
    my $mask = 2**$jj;
    if (($blade & $mask)) {
      $runtot++;
    }
  }
  $map[$dim] = $runtot;
  return @map;
}

sub clifmul{
  my ($A, $B, $dotwedge) = @_;
  if (! defined $dotwedge){
    $dotwedge = 0;
  }
  if (ref($A) ne 'CLIF' || ref($B) ne 'CLIF') {
    die "clifmul only works on new CLIF representation.\n";
  }
  my $Acomp = @$A;
  my $Bcomp = @$B;
  my $ABcomp = $Acomp;
  if ($Bcomp > $Acomp){
    $ABcomp = $Bcomp;             # whichever is greater
  }
  my $rslt = [];                # new anonymous array
  for (my $ii = 0; $ii < $ABcomp; $ii++) {
    @$rslt[$ii] = 0;
  }
# Note that $ii and $jj are bitcodes
  for (my $ii = 0; $ii < $Acomp; $ii++) {
    for (my $jj = 0; $jj < $Bcomp; $jj++) {
# inner loop over all blade*blade products
      my $Agrade = bladegrade($ii);
      my $Bgrade = bladegrade($jj);
      my $goodgrade = -9e99;
      if ($dotwedge == 1) {
        $goodgrade = $Agrade + $Bgrade;
      } elsif ($dotwedge == -1) {
        $goodgrade = abs($Agrade - $Bgrade);
      }

      my $prod = $ii ^ $jj;
      if ($goodgrade < 0 || bladegrade($prod) == $goodgrade) {
        @$rslt[$prod] += sign_of_product($ii, $jj)
          * @$A[$ii] * @$B[$jj];
      }
    }
  }
  bless $rslt, 'CLIF';
  return $rslt;
}

sub promote{
  my ($A) = @_;
  my $clif;
  my $type = myref($A);
  if ($type eq 'CLIF') {
    return $A;
  }
  if ($type eq 'plain_number') {
# clif in a zero-dimensional space; no basis vectors:
    $clif = [$A];
  }

# New D=3 clif packing rule
# 0, 1, 2, 3,  4, 5!,  6,  7
# .  x, y, xy, z, xz, yz, xyz
  elsif ($type eq 'QUAT') {
    my ($As, $Ayz, $Azx, $Axy) = @$A;
    $clif = [$As, 0, 0, $Axy,
                 0, -$Azx, $Ayz, 0];
  }
  elsif ($type eq 'VECTOR') {
    my $dim = @$A;
    my $clifsize = 2**$dim;
    for (my $ii = 0; $ii < $clifsize; $ii++){
      @$clif[$ii] = 0;
    }
    my $mask = 1;
    foreach my $comp (@$A) {
      @$clif[$mask] = $comp;
      $mask <<= 1;
    }
  }
  else {
    die "Don't know how to promote object of type '$type'.\n";
  }
  bless $clif, 'CLIF';
  return $clif;
}

sub add{
  my ($A, $B) = @_;
  if (ref($A) ne ref($B)){      # mismatch; promote them both
    $A = promote($A);
    $B = promote($B);
  }
  if (myref($A) eq 'plain_number') {
    return $A + $B;
  }
  my $rslt;
  my $smallguy;
  if (@$A > @$B) {
    $rslt = copy($A);
    $smallguy = $B;
  } else {
    $rslt = copy($B);
    $smallguy = $A;
  }
  my $ncomp = @$smallguy;       # the lesser of two weevils
  for (my $ii = 0; $ii < $ncomp; $ii++){
    @$rslt[$ii] += @$smallguy[$ii];
  }
  bless $rslt, ref($A);
  return $rslt;
}

sub copy{
  my ($A) = @_;
  my $ncomp = @$A;              # number of components
  my $C = [];                   # new anonymous array
  for (my $ii = 0; $ii < $ncomp; $ii++){
    @$C[$ii] = @$A[$ii];
  }
  bless $C, ref($A);
  return $C;
}

1;