marp.pl

use strict;
use warnings;

use 5.010;
use PDL;    # must be called before (!) 'use Inline Pdlpp' calls
use IO::All -binary;
use Inline 'Pdlpp' => io( "pdlpp.pm" )->all;
use Benchmark 'timethis';

my $v = pdl( 1, 0, 0 );
$v->rotate( 180, pdl( 0, 1, 0 ) );
say $v;

my $vec  = pdl( 1, 0, 0 );
my $axis = pdl( 0, 1, 0 );
timethis( 40000, sub { $vec->rotate( 180, $axis ) } );

pdlpp.pm

use strict;


pp_def(
    'rotate',
    Pars => 'vec(n); angle(); axis(n)',
    Code => q{
        double PI = 3.141592653589793;
        double SinHalfAngle = sin( $angle() / 2.0 / 180.0 * PI );
        double CosHalfAngle = cos( $angle() / 2.0 / 180.0 * PI );
        
        // rotation quaternion
        double Rx = $axis(n=>0) * SinHalfAngle;
        double Ry = $axis(n=>1) * SinHalfAngle;
        double Rz = $axis(n=>2) * SinHalfAngle;
        double Rw = CosHalfAngle;
        
        double Vx = $vec(n=>0);
        double Vy = $vec(n=>1);
        double Vz = $vec(n=>2);
        double Vw = 0;
        
        // first step
        double x1 = ( Rx * Vw ) + ( Rw * Vx ) + ( Ry * Vz ) - ( Rz * Vy );
        double y1 = ( Ry * Vw ) + ( Rw * Vy ) + ( Rz * Vx ) - ( Rx * Vz );
        double z1 = ( Rz * Vw ) + ( Rw * Vz ) + ( Rx * Vy ) - ( Ry * Vx );
        double w1 = ( Rw * Vw ) - ( Rx * Vx ) - ( Ry * Vy ) - ( Rz * Vz );

        // conjugated quaternion
        double Qx = Rx * -1;
        double Qy = Ry * -1;
        double Qz = Rz * -1;
        double Qw = CosHalfAngle;
        
        // second step
        double x2 = ( x1 * Qw ) + ( w1 * Qx ) + ( y1 * Qz ) - ( z1 * Qy );
        double y2 = ( y1 * Qw ) + ( w1 * Qy ) + ( z1 * Qx ) - ( x1 * Qz );
        double z2 = ( z1 * Qw ) + ( w1 * Qz ) + ( x1 * Qy ) - ( y1 * Qx );
        double w2 = ( w1 * Qw ) - ( x1 * Qx ) - ( y1 * Qy ) - ( z1 * Qz );
        
        $vec(n=>0) = x2;
        $vec(n=>1) = y2;
        $vec(n=>2) = z2;
    },
);