home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 505a.lha / GrapicsGems / AALines / FastMatMul.c < prev    next >
Text File  |  1991-05-01  |  12KB  |  509 lines
  1. /*  FILENAME: FastMatMul.c  [revised 18 AUG 90]
  2.  
  3.     AUTHOR:  Kelvin Thompson
  4.  
  5.     DESCRIPTION:  Routines to multiply different kinds of 4x4
  6.       matrices as fast as possible.  Based on ideas on pages 454,
  7.       460-461, and 646 of _Graphics_Gems_.
  8.  
  9.     These routines offer two advantages over the standard
  10.       V3MatMul in the _Graphics_Gems_ vector library GGVecLib.c.
  11.       First, the routines are faster.  Second, they can handle
  12.       input matrices that are the same as the output matrix.
  13.       The routines have the disadvantage of taking more code
  14.       space (from unwound loops).
  15.  
  16.         The routines are about as fast as you can get for
  17.       general-purpose multiplication.  If you have special
  18.       knowledge about your system, you may be able to improve
  19.       them a little more:
  20.  
  21.         [1]  If you know that your input and output matrices will
  22.       never overlap: remove the tests near the beginning and end of
  23.       each routine, and just #define 'mptr' to 'result'.  (The
  24.       standard library's V3MatMul makes this assumption.)
  25.  
  26.     [2]  If you know that your compiler supports more than
  27.       three pointer-to-double registers in a subroutine: make
  28.       'result' in each routine a register variable.  You might
  29.       also make the 'usetemp' boolean a register.
  30.  
  31.     [3]  If you have limited stack space, or your system can
  32.       access global memory faster than stack:  make each routine's
  33.       'tempx' a static, or let all routines share a global 'tempx'.
  34.       (This is useless if assumption [1] holds.)
  35. */
  36.  
  37. /* definitions from "GraphicsGems.h" */
  38. typedef struct Matrix3Struct {    /* 3-by-3 matrix */
  39.     double element[3][3];
  40.     } Matrix3;
  41. typedef struct Matrix4Struct {    /* 4-by-4 matrix */
  42.     double element[4][4];
  43.     } Matrix4;
  44.  
  45. /* routines in this file */
  46. Matrix3 *V2MatMul();     /* general 3x3 matrix multiplier */
  47. Matrix4 *V3MatMul();     /* general 4x4 matrix multiplier */
  48. Matrix4 *V3AffMatMul();  /* affine 4x4 matrix multiplier */
  49. Matrix4 *V3LinMatMul();  /* linear 4x4 matrix multiplier */
  50.  
  51. /* macro to access matrix element */
  52. #define MVAL(mptr,row,col)  ((mptr)->element[row][col])
  53.  
  54.  
  55.  
  56.  
  57.  
  58. /*  ************************************
  59.     *                                  *
  60.     *           V2MatMul               *
  61.     *                                  *
  62.     ************************************
  63.  
  64.     DESCRIPTION:  Multiply two general 3x3 matricies.  If one of
  65.       the input matrices is the same as the output, write the
  66.       result to a temporary matrix during multiplication, then
  67.       copy to the output matrix.
  68.  
  69.     ENTRY:
  70.       a -- pointer to left matrix
  71.       b -- pointer to right matrix
  72.       result -- result matrix
  73.  
  74.     EXIT:  returns 'result'
  75. */
  76.  
  77. Matrix3 *V2MatMul ( a, b, result )
  78. register Matrix3 *a,*b;
  79. Matrix3 *result;
  80. {
  81. register Matrix3 *mptr;
  82. int usetemp;  /* boolean */
  83. Matrix3 tempx;
  84.  
  85. /* decide where intermediate result goes */
  86. usetemp = ( a == result  ||  b == result );
  87. if ( usetemp )
  88.   mptr = & tempx;
  89. else
  90.   mptr = result;
  91.  
  92. MVAL(mptr,0,0) =
  93.      MVAL(a,0,0) * MVAL(b,0,0)
  94.   +  MVAL(a,0,1) * MVAL(b,1,0)
  95.   +  MVAL(a,0,2) * MVAL(b,2,0);
  96.  
  97. MVAL(mptr,0,1) =
  98.      MVAL(a,0,0) * MVAL(b,0,1)
  99.   +  MVAL(a,0,1) * MVAL(b,1,1)
  100.   +  MVAL(a,0,2) * MVAL(b,2,1);
  101.  
  102. MVAL(mptr,0,2) =
  103.      MVAL(a,0,0) * MVAL(b,0,2)
  104.   +  MVAL(a,0,1) * MVAL(b,1,2)
  105.   +  MVAL(a,0,2) * MVAL(b,2,2);
  106.  
  107. MVAL(mptr,1,0) =
  108.      MVAL(a,1,0) * MVAL(b,0,0)
  109.   +  MVAL(a,1,1) * MVAL(b,1,0)
  110.   +  MVAL(a,1,2) * MVAL(b,2,0);
  111.  
  112. MVAL(mptr,1,1) =
  113.      MVAL(a,1,0) * MVAL(b,0,1)
  114.   +  MVAL(a,1,1) * MVAL(b,1,1)
  115.   +  MVAL(a,1,2) * MVAL(b,2,1);
  116.  
  117. MVAL(mptr,1,2) =
  118.      MVAL(a,1,0) * MVAL(b,0,2)
  119.   +  MVAL(a,1,1) * MVAL(b,1,2)
  120.   +  MVAL(a,1,2) * MVAL(b,2,2);
  121.  
  122. MVAL(mptr,2,0) =
  123.      MVAL(a,2,0) * MVAL(b,0,0)
  124.   +  MVAL(a,2,1) * MVAL(b,1,0)
  125.   +  MVAL(a,2,2) * MVAL(b,2,0);
  126.  
  127. MVAL(mptr,2,1) =
  128.      MVAL(a,2,0) * MVAL(b,0,1)
  129.   +  MVAL(a,2,1) * MVAL(b,1,1)
  130.   +  MVAL(a,2,2) * MVAL(b,2,1);
  131.  
  132. MVAL(mptr,2,2) =
  133.      MVAL(a,2,0) * MVAL(b,0,2)
  134.   +  MVAL(a,2,1) * MVAL(b,1,2)
  135.   +  MVAL(a,2,2) * MVAL(b,2,2);
  136.  
  137. /* copy temp matrix to result if needed */
  138. if ( usetemp )
  139.   *result = *mptr;
  140.  
  141. return result;
  142. }
  143.  
  144.  
  145.  
  146.  
  147. /*  ************************************
  148.     *                                  *
  149.     *           V3MatMul               *
  150.     *                                  *
  151.     ************************************
  152.  
  153.     DESCRIPTION:  Multiply two general 4x4 matricies.  If one of
  154.       the input matrices is the same as the output, write the
  155.       result to a temporary matrix during multiplication, then
  156.       copy to the output matrix.
  157.  
  158.     ENTRY:
  159.       a -- pointer to left matrix
  160.       b -- pointer to right matrix
  161.       result -- result matrix
  162.  
  163.     EXIT:  returns 'result'
  164. */
  165.  
  166. Matrix4 *V3MatMul ( a, b, result )
  167. register Matrix4 *a,*b;
  168. Matrix4 *result;
  169. {
  170. register Matrix4 *mptr;
  171. int usetemp;  /* boolean */
  172. Matrix4 tempx;
  173.  
  174. /* decide where intermediate result goes */
  175. usetemp = ( a == result  ||  b == result );
  176. if ( usetemp )
  177.   mptr = & tempx;
  178. else
  179.   mptr = result;
  180.  
  181. MVAL(mptr,0,0) =
  182.      MVAL(a,0,0) * MVAL(b,0,0)
  183.   +  MVAL(a,0,1) * MVAL(b,1,0)
  184.   +  MVAL(a,0,2) * MVAL(b,2,0)
  185.   +  MVAL(a,0,3) * MVAL(b,3,0);
  186.  
  187. MVAL(mptr,0,1) =
  188.      MVAL(a,0,0) * MVAL(b,0,1)
  189.   +  MVAL(a,0,1) * MVAL(b,1,1)
  190.   +  MVAL(a,0,2) * MVAL(b,2,1)
  191.   +  MVAL(a,0,3) * MVAL(b,3,1);
  192.  
  193. MVAL(mptr,0,2) =
  194.      MVAL(a,0,0) * MVAL(b,0,2)
  195.   +  MVAL(a,0,1) * MVAL(b,1,2)
  196.   +  MVAL(a,0,2) * MVAL(b,2,2)
  197.   +  MVAL(a,0,3) * MVAL(b,3,2);
  198.  
  199. MVAL(mptr,0,3) =
  200.      MVAL(a,0,0) * MVAL(b,0,3)
  201.   +  MVAL(a,0,1) * MVAL(b,1,3)
  202.   +  MVAL(a,0,2) * MVAL(b,2,3)
  203.   +  MVAL(a,0,3) * MVAL(b,3,3);
  204.  
  205. MVAL(mptr,1,0) =
  206.      MVAL(a,1,0) * MVAL(b,0,0)
  207.   +  MVAL(a,1,1) * MVAL(b,1,0)
  208.   +  MVAL(a,1,2) * MVAL(b,2,0)
  209.   +  MVAL(a,1,3) * MVAL(b,3,0);
  210.  
  211. MVAL(mptr,1,1) =
  212.      MVAL(a,1,0) * MVAL(b,0,1)
  213.   +  MVAL(a,1,1) * MVAL(b,1,1)
  214.   +  MVAL(a,1,2) * MVAL(b,2,1)
  215.   +  MVAL(a,1,3) * MVAL(b,3,1);
  216.  
  217. MVAL(mptr,1,2) =
  218.      MVAL(a,1,0) * MVAL(b,0,2)
  219.   +  MVAL(a,1,1) * MVAL(b,1,2)
  220.   +  MVAL(a,1,2) * MVAL(b,2,2)
  221.   +  MVAL(a,1,3) * MVAL(b,3,2);
  222.  
  223. MVAL(mptr,1,3) =
  224.      MVAL(a,1,0) * MVAL(b,0,3)
  225.   +  MVAL(a,1,1) * MVAL(b,1,3)
  226.   +  MVAL(a,1,2) * MVAL(b,2,3)
  227.   +  MVAL(a,1,3) * MVAL(b,3,3);
  228.  
  229. MVAL(mptr,2,0) =
  230.      MVAL(a,2,0) * MVAL(b,0,0)
  231.   +  MVAL(a,2,1) * MVAL(b,1,0)
  232.   +  MVAL(a,2,2) * MVAL(b,2,0)
  233.   +  MVAL(a,2,3) * MVAL(b,3,0);
  234.  
  235. MVAL(mptr,2,1) =
  236.      MVAL(a,2,0) * MVAL(b,0,1)
  237.   +  MVAL(a,2,1) * MVAL(b,1,1)
  238.   +  MVAL(a,2,2) * MVAL(b,2,1)
  239.   +  MVAL(a,2,3) * MVAL(b,3,1);
  240.  
  241. MVAL(mptr,2,2) =
  242.      MVAL(a,2,0) * MVAL(b,0,2)
  243.   +  MVAL(a,2,1) * MVAL(b,1,2)
  244.   +  MVAL(a,2,2) * MVAL(b,2,2)
  245.   +  MVAL(a,2,3) * MVAL(b,3,2);
  246.  
  247. MVAL(mptr,2,3) =
  248.      MVAL(a,2,0) * MVAL(b,0,3)
  249.   +  MVAL(a,2,1) * MVAL(b,1,3)
  250.   +  MVAL(a,2,2) * MVAL(b,2,3)
  251.   +  MVAL(a,2,3) * MVAL(b,3,3);
  252.  
  253. MVAL(mptr,3,0) =
  254.      MVAL(a,3,0) * MVAL(b,0,0)
  255.   +  MVAL(a,3,1) * MVAL(b,1,0)
  256.   +  MVAL(a,3,2) * MVAL(b,2,0)
  257.   +  MVAL(a,3,3) * MVAL(b,3,0);
  258.  
  259. MVAL(mptr,3,1) =
  260.      MVAL(a,3,0) * MVAL(b,0,1)
  261.   +  MVAL(a,3,1) * MVAL(b,1,1)
  262.   +  MVAL(a,3,2) * MVAL(b,2,1)
  263.   +  MVAL(a,3,3) * MVAL(b,3,1);
  264.  
  265. MVAL(mptr,3,2) =
  266.      MVAL(a,3,0) * MVAL(b,0,2)
  267.   +  MVAL(a,3,1) * MVAL(b,1,2)
  268.   +  MVAL(a,3,2) * MVAL(b,2,2)
  269.   +  MVAL(a,3,3) * MVAL(b,3,2);
  270.  
  271. MVAL(mptr,3,3) =
  272.      MVAL(a,3,0) * MVAL(b,0,3)
  273.   +  MVAL(a,3,1) * MVAL(b,1,3)
  274.   +  MVAL(a,3,2) * MVAL(b,2,3)
  275.   +  MVAL(a,3,3) * MVAL(b,3,3);
  276.  
  277. /* copy temp matrix to result if needed */
  278. if ( usetemp )
  279.   *result = *mptr;
  280.  
  281. return result;
  282. }
  283.  
  284.  
  285.  
  286.  
  287.  
  288. /*  ************************************
  289.     *                                  *
  290.     *        V3AffMatMul               *
  291.     *                                  *
  292.     ************************************
  293.  
  294.     DESCRIPTION:  Multiply two affine 4x4 matricies.  The
  295.       routine assumes the rightmost column of each input
  296.       matrix is [0 0 0 1].  The output matrix will have the
  297.       same property.
  298.     
  299.       If one of the input matrices is the same as the output,
  300.       write the result to a temporary matrix during multiplication,
  301.       then copy to the output matrix.
  302.  
  303.     ENTRY:
  304.       a -- pointer to left matrix
  305.       b -- pointer to right matrix
  306.       result -- result matrix
  307.  
  308.     EXIT:  returns 'result'
  309. */
  310.  
  311. Matrix4 *V3AffMatMul ( a, b, result )
  312. register Matrix4 *a,*b;
  313. Matrix4 *result;
  314. {
  315. register Matrix4 *mptr;
  316. int usetemp;  /* boolean */
  317. Matrix4 tempx;
  318.  
  319. /* decide where intermediate result goes */
  320. usetemp = ( a == result  ||  b == result );
  321. if ( usetemp )
  322.   mptr = & tempx;
  323. else
  324.   mptr = result;
  325.  
  326. MVAL(mptr,0,0) =
  327.      MVAL(a,0,0) * MVAL(b,0,0)
  328.   +  MVAL(a,0,1) * MVAL(b,1,0)
  329.   +  MVAL(a,0,2) * MVAL(b,2,0);
  330.  
  331. MVAL(mptr,0,1) =
  332.      MVAL(a,0,0) * MVAL(b,0,1)
  333.   +  MVAL(a,0,1) * MVAL(b,1,1)
  334.   +  MVAL(a,0,2) * MVAL(b,2,1);
  335.  
  336. MVAL(mptr,0,2) =
  337.      MVAL(a,0,0) * MVAL(b,0,2)
  338.   +  MVAL(a,0,1) * MVAL(b,1,2)
  339.   +  MVAL(a,0,2) * MVAL(b,2,2);
  340.  
  341. MVAL(mptr,0,3) = 0.0;
  342.  
  343. MVAL(mptr,1,0) =
  344.      MVAL(a,1,0) * MVAL(b,0,0)
  345.   +  MVAL(a,1,1) * MVAL(b,1,0)
  346.   +  MVAL(a,1,2) * MVAL(b,2,0);
  347.  
  348. MVAL(mptr,1,1) =
  349.      MVAL(a,1,0) * MVAL(b,0,1)
  350.   +  MVAL(a,1,1) * MVAL(b,1,1)
  351.   +  MVAL(a,1,2) * MVAL(b,2,1);
  352.  
  353. MVAL(mptr,1,2) =
  354.      MVAL(a,1,0) * MVAL(b,0,2)
  355.   +  MVAL(a,1,1) * MVAL(b,1,2)
  356.   +  MVAL(a,1,2) * MVAL(b,2,2);
  357.  
  358. MVAL(mptr,1,3) = 0.0;
  359.  
  360. MVAL(mptr,2,0) =
  361.      MVAL(a,2,0) * MVAL(b,0,0)
  362.   +  MVAL(a,2,1) * MVAL(b,1,0)
  363.   +  MVAL(a,2,2) * MVAL(b,2,0);
  364.  
  365. MVAL(mptr,2,1) =
  366.      MVAL(a,2,0) * MVAL(b,0,1)
  367.   +  MVAL(a,2,1) * MVAL(b,1,1)
  368.   +  MVAL(a,2,2) * MVAL(b,2,1);
  369.  
  370. MVAL(mptr,2,2) =
  371.      MVAL(a,2,0) * MVAL(b,0,2)
  372.   +  MVAL(a,2,1) * MVAL(b,1,2)
  373.   +  MVAL(a,2,2) * MVAL(b,2,2);
  374.  
  375. MVAL(mptr,2,3) = 0.0;
  376.  
  377. MVAL(mptr,3,0) =
  378.      MVAL(a,3,0) * MVAL(b,0,0)
  379.   +  MVAL(a,3,1) * MVAL(b,1,0)
  380.   +  MVAL(a,3,2) * MVAL(b,2,0)
  381.   +                MVAL(b,3,0);
  382.  
  383. MVAL(mptr,3,1) =
  384.      MVAL(a,3,0) * MVAL(b,0,1)
  385.   +  MVAL(a,3,1) * MVAL(b,1,1)
  386.   +  MVAL(a,3,2) * MVAL(b,2,1)
  387.   +                MVAL(b,3,1);
  388.  
  389. MVAL(mptr,3,2) =
  390.      MVAL(a,3,0) * MVAL(b,0,2)
  391.   +  MVAL(a,3,1) * MVAL(b,1,2)
  392.   +  MVAL(a,3,2) * MVAL(b,2,2)
  393.   +                MVAL(b,3,2);
  394.  
  395. MVAL(mptr,3,3) = 1.0;
  396.  
  397. /* copy temp matrix to result if needed */
  398. if ( usetemp )
  399.   *result = *mptr;
  400.  
  401. return result;
  402. }
  403.  
  404.  
  405.  
  406.  
  407.  
  408. /*  ************************************
  409.     *                                  *
  410.     *        V3LinMatMul               *
  411.     *                                  *
  412.     ************************************
  413.  
  414.     DESCRIPTION:  Multiply two affine 4x4 matricies.  The
  415.       routine assumes the right column and bottom line
  416.       of each input matrix is [0 0 0 1].  The output matrix
  417.       will have the same property.  This is pretty much the
  418.       same thing as multiplying two 3x3 matrices.
  419.     
  420.       If one of the input matrices is the same as the output,
  421.       write the result to a temporary matrix during multiplication,
  422.       then copy to the output matrix.
  423.  
  424.     ENTRY:
  425.       a -- pointer to left matrix
  426.       b -- pointer to right matrix
  427.       result -- result matrix
  428.  
  429.     EXIT:  returns 'result'
  430. */
  431.  
  432. Matrix4 *V3LinMatMul ( a, b, result )
  433. register Matrix4 *a,*b;
  434. Matrix4 *result;
  435. {
  436. register Matrix4 *mptr;
  437. int usetemp;  /* boolean */
  438. Matrix4 tempx;
  439.  
  440. /* decide where intermediate result goes */
  441. usetemp = ( a == result  ||  b == result );
  442. if ( usetemp )
  443.   mptr = & tempx;
  444. else
  445.   mptr = result;
  446.  
  447. MVAL(mptr,0,0) =
  448.      MVAL(a,0,0) * MVAL(b,0,0)
  449.   +  MVAL(a,0,1) * MVAL(b,1,0)
  450.   +  MVAL(a,0,2) * MVAL(b,2,0);
  451.  
  452. MVAL(mptr,0,1) =
  453.      MVAL(a,0,0) * MVAL(b,0,1)
  454.   +  MVAL(a,0,1) * MVAL(b,1,1)
  455.   +  MVAL(a,0,2) * MVAL(b,2,1);
  456.  
  457. MVAL(mptr,0,2) =
  458.      MVAL(a,0,0) * MVAL(b,0,2)
  459.   +  MVAL(a,0,1) * MVAL(b,1,2)
  460.   +  MVAL(a,0,2) * MVAL(b,2,2);
  461.  
  462. MVAL(mptr,0,3) = 0.0;
  463.  
  464. MVAL(mptr,1,0) =
  465.      MVAL(a,1,0) * MVAL(b,0,0)
  466.   +  MVAL(a,1,1) * MVAL(b,1,0)
  467.   +  MVAL(a,1,2) * MVAL(b,2,0);
  468.  
  469. MVAL(mptr,1,1) =
  470.      MVAL(a,1,0) * MVAL(b,0,1)
  471.   +  MVAL(a,1,1) * MVAL(b,1,1)
  472.   +  MVAL(a,1,2) * MVAL(b,2,1);
  473.  
  474. MVAL(mptr,1,2) =
  475.      MVAL(a,1,0) * MVAL(b,0,2)
  476.   +  MVAL(a,1,1) * MVAL(b,1,2)
  477.   +  MVAL(a,1,2) * MVAL(b,2,2);
  478.  
  479. MVAL(mptr,1,3) = 0.0;
  480.  
  481. MVAL(mptr,2,0) =
  482.      MVAL(a,2,0) * MVAL(b,0,0)
  483.   +  MVAL(a,2,1) * MVAL(b,1,0)
  484.   +  MVAL(a,2,2) * MVAL(b,2,0);
  485.  
  486. MVAL(mptr,2,1) =
  487.      MVAL(a,2,0) * MVAL(b,0,1)
  488.   +  MVAL(a,2,1) * MVAL(b,1,1)
  489.   +  MVAL(a,2,2) * MVAL(b,2,1);
  490.  
  491. MVAL(mptr,2,2) =
  492.      MVAL(a,2,0) * MVAL(b,0,2)
  493.   +  MVAL(a,2,1) * MVAL(b,1,2)
  494.   +  MVAL(a,2,2) * MVAL(b,2,2);
  495.  
  496. MVAL(mptr,2,3) = 0.0;
  497.  
  498. MVAL(mptr,3,0) = 0.0;
  499. MVAL(mptr,3,1) = 0.0;
  500. MVAL(mptr,3,2) = 0.0;
  501. MVAL(mptr,3,3) = 1.0;
  502.  
  503. /* copy temp matrix to result if needed */
  504. if ( usetemp )
  505.   *result = *mptr;
  506.  
  507. return result;
  508. }
  509.