Source code in the pl/i programming language

/* dot product, cross product, etc.                  4 June 2011 */

test_products: procedure options (main);

   declare a(3) fixed initial (3, 4, 5);
   declare b(3) fixed initial (4, 3, 5);
   declare c(3) fixed initial (-5, -12, -13);
   declare e(3) fixed;

   put skip list ('a . b =', dot_product(a, b));
   call cross_product(a, b, e);  put skip list ('a x b =', e);
   put skip list ('a . (b x c) =',  scalar_triple_product(a, b, c));
   call vector_triple_product(a, b, c, e); put skip list ('a x (b x c) =', e);


dot_product: procedure (a, b) returns (fixed);
   declare (a, b) (*) fixed;
   return (sum(a*b));
end dot_product;

cross_product: procedure (a, b, c);
   declare (a, b, c) (*) fixed;
   c(1) = a(2)*b(3) - a(3)*b(2);
   c(2) = a(3)*b(1) - a(1)*b(3);
   c(3) = a(1)*b(2) - a(2)*b(1);
end cross_product;

scalar_triple_product: procedure (a, b, c) returns (fixed);
   declare (a, b, c)(*) fixed;
   declare t(hbound(a, 1)) fixed;
   call cross_product(b, c, t);
   return (dot_product(a, t));
end scalar_triple_product;

vector_triple_product: procedure (a, b, c, e);
   declare (a, b, c, e)(*) fixed;
   declare t(hbound(a,1))  fixed;
   call cross_product(b, c, t);
   call cross_product(a, t, e);
end vector_triple_product;

end test_products;

/* This version uses the ability of PL/I to return arrays. */

/* dot product, cross product, etc.                  6 June 2011 */

test_products: procedure options (main);
   define structure 1 vector, 2 vec(3) fixed;
   declare (a, b, c) type(vector);

   a.vec(1) =  3; a.vec(2) =   4; a.vec(3) =   5;
   b.vec(1) =  4; b.vec(2) =   3; b.vec(3) =   5;
   c.vec(1) = -5; c.vec(2) = -12; c.vec(3) = -13;

   put skip list ('a . b =',       dot_product  (a, b) );
   put skip list ('a x b =',       cross_product(a, b).vec);
   put skip list ('a . (b x c) =', scalar_triple_product(a, b, c) );
   put skip list ('a x (b x c) =', vector_triple_product(a, b, c).vec);


dot_product: procedure (a, b) returns (fixed);
   declare (a, b) type(vector);
   return (sum(a.vec*b.vec));
end dot_product;

cross_product: procedure (a, b) returns (type(vector));
   declare (a, b) type(vector);
   declare c type vector;
   c.vec(1) = a.vec(2)*b.vec(3) - a.vec(3)*b.vec(2);
   c.vec(2) = a.vec(3)*b.vec(1) - a.vec(1)*b.vec(3);
   c.vec(3) = a.vec(1)*b.vec(2) - a.vec(2)*b.vec(1);
   return (c);
end cross_product;

scalar_triple_product: procedure (a, b, c) returns (fixed);
   declare (a, b, c) type(vector);
   declare t type (vector);
   t =  cross_product(b, c);
   return (dot_product(a, t));
end scalar_triple_product;

vector_triple_product: procedure (a, b, c) returns (type(vector));
   declare (a, b, c) type(vector);
   declare (t, e) type (vector);
   t = cross_product(b, c);
   e = cross_product(a, t);
   return (e);
end vector_triple_product;

end test_products;