Trans.cpp
上传用户:lhwx1029
上传日期:2013-03-07
资源大小:1173k
文件大小:20k
源码类别:

3D图形编程

开发平台:

Visual C++

  1. /** 3DGPL *************************************************
  2.  * ()                                                     *
  3.  * 3-D transformations of coordinates.                    *
  4.  *                                                        *
  5.  * Ifdefs:                                                *
  6.  *  _FLOAT_                  Use float implementation;    *
  7.  *  _FIXED_                  Use fixed implementation;    *
  8.  *  _Z_BUFFER_               Depth array;                 *
  9.  *  _PAINTER_                Back to front order.         *
  10.  *                                                        *
  11.  * Defines:                                               *
  12.  *  T_init_math              Creating sin/cos tables;     *
  13.  *                                                        *
  14.  *  T_translation            Translation of coordinates;  *
  15.  *  T_scaling                Scaling coordinates;         *
  16.  *                                                        *
  17.  *  T_set_self_rotation      Setting object rotation;     *
  18.  *  T_self_rotation          Rotation of coordinates;     *
  19.  *  T_set_world_rotation     Setting viewer's rotation;   *
  20.  *  T_world_rotation         Rotation of coordinates;     *
  21.  *  T_cancatinate_self_world The two above matrices;      *
  22.  *  T_cancatinated_rotation  From object into view space; *
  23.  *                                                        *
  24.  *  T_perspective            Transform to perspective;    *
  25.  *                                                        *
  26.  *  T_linear_solve           Gaussian elimination.        *
  27.  *                                                        *
  28.  * (c) 1995-98 Sergei Savchenko, (savs@cs.mcgill.ca)      *
  29. **********************************************************/
  30. #include "RayTracing.h"           /* fast copy routines */
  31. #include "Trans.h"                 /* 3D mathematics */
  32. #include <math.h>                           /* sin & cos */
  33. #define T_RADS 40.7436642                   /* pseudo-grads into rads */
  34. int T_log_focus;                            /* foreshortening */
  35. #if defined(_FIXED_)
  36. typedef HW_32_bit T_fixed;                  /* better be 32bit machine */
  37. T_fixed T_wx1,T_wx2,T_wx3,T_wy1,T_wy2,T_wy3,T_wz1,T_wz2,T_wz3;
  38. T_fixed T_sx1,T_sx2,T_sx3,T_sy1,T_sy2,T_sy3,T_sz1,T_sz2,T_sz3;
  39. T_fixed T_cx1,T_cx2,T_cx3,T_cy1,T_cy2,T_cy3,T_cz1,T_cz2,T_cz3;
  40. T_fixed T_tx,T_ty,T_tz;
  41. T_fixed T_sin[256],T_cos[256];              /* precalculated */
  42. #endif
  43. #if defined(_FLOAT_)
  44. float T_wx1,T_wx2,T_wx3,T_wy1,T_wy2,T_wy3,T_wz1,T_wz2,T_wz3;
  45. float T_sx1,T_sx2,T_sx3,T_sy1,T_sy2,T_sy3,T_sz1,T_sz2,T_sz3;
  46. float T_cx1,T_cx2,T_cx3,T_cy1,T_cy2,T_cy3,T_cz1,T_cz2,T_cz3;
  47. float T_tx,T_ty,T_tz;
  48. float T_sin[256],T_cos[256];                /* precalculated */
  49. #endif
  50. /**********************************************************
  51.  * Initializing tables of trigonometric functions.        *
  52.  *                                                        *
  53.  * SETS: T_sin and T_cos arrays.                          *
  54.  * -----                                                  *
  55. **********************************************************/
  56. void T_init_math(void)
  57. {
  58.  int i;
  59.  for(i=0;i<256;i++)
  60.  {
  61. #if defined(_FIXED_)
  62.   T_sin[i]=sin(i/T_RADS)*(1<<T_P);          /* moving fraction into int part */
  63.   T_cos[i]=cos(i/T_RADS)*(1<<T_P);          /* which is cast into fixed */
  64. #endif
  65. #if defined(_FLOAT_)
  66.   T_sin[i]=sin(i/T_RADS);
  67.   T_cos[i]=cos(i/T_RADS);
  68. #endif
  69.  }
  70. }
  71. /**********************************************************
  72.  * Translation of coordinates.                            *
  73. **********************************************************/
  74. //register: the variable is to be stored in a machine register, if possible
  75. void T_translation(int *from,register int *to,int length,
  76.                    int addx,int addy,int addz
  77.                   )
  78. {
  79.  register int i;
  80.  for(i=0;i<length;i++)
  81.  {
  82.   *to++=(*from++)+addx;
  83.   *to++=(*from++)+addy;
  84.   *to++=(*from++)+addz;                     /* translating X Y Z */
  85.   to++; from++;                             /* 4th element */
  86. //the 4th element saves for what?
  87.  }
  88. }
  89. /**********************************************************
  90.  * Scaling of coordinates.                                *
  91. **********************************************************/
  92. void T_scaling(int *from,register int *to,int length,
  93.                int mulx,int muly,int mulz
  94.               )
  95. {
  96.  register int i;
  97.  for(i=0;i<length;i++)
  98.  {
  99.   *to++=(*from++)*mulx;
  100.   *to++=(*from++)*muly;
  101.   *to++=(*from++)*mulz;                     /* scaling X Y Z */
  102.   to++; from++;                             /* 4th element */
  103.  }
  104. }
  105. /**********************************************************
  106.  * Constructing rotation matrix for object to world       *
  107.  * rotation: (alp-bet-gam), roll then pitch then yaw      *
  108.  * sequence. gam-yaw, bet-pitch, alp-roll. NB! angles     *
  109.  * assumed to be in pseudo-degrees in the range [0,256).  *
  110.  *                                                        *
  111.  *          Y^    Z           x'=y*sin(alp)+x*cos(alp)    *
  112.  *           |   /            y'=y*cos(alp)-x*sin(alp)    *
  113.  *           |  / alp         z'=z                        *
  114.  *          /|<---+                                       *
  115.  *         | |/   |           x"=x'                       *
  116.  *  -------|-+-----------> X  y"=y'*cos(bet)-z'*sin(bet)  *
  117.  *     bet | |   __           z"=y'*sin(bet)+z'*cos(bet)  *
  118.  *         V/|   /| gam                                   *
  119.  *         /----+             x"'=z"*sin(gam)+x"*cos(gam) *
  120.  *        /  |                y"'=y"                      *
  121.  *       /   |                z"'=z"*cos(gam)-x"*sin(gam) *
  122.  *           |                                            *
  123.  *                                                        *
  124.  * SETS: T_sx1,...,T_sz3.                                 *
  125.  * -----                                                  *
  126. **********************************************************/
  127. void T_set_self_rotation(unsigned char alp,unsigned char bet,
  128.                          unsigned char gam)
  129. {
  130. #if defined(_FIXED_)
  131.  T_fixed cosalp,sinalp,cosbet,sinbet,cosgam,singam;
  132. #endif
  133. #if defined(_FLOAT_)
  134.  float cosalp,sinalp,cosbet,sinbet,cosgam,singam;
  135. #endif
  136.  cosalp=T_cos[alp];
  137.  sinalp=T_sin[alp];
  138.  cosbet=T_cos[bet];
  139.  sinbet=T_sin[bet];
  140.  cosgam=T_cos[gam];
  141.  singam=T_sin[gam];                         /* initializing */
  142. #if defined(_FIXED_)
  143.  T_sx1=((cosalp*cosgam)>>T_P)-((sinalp*((sinbet*singam)>>T_P))>>T_P);
  144.  T_sy1=((sinalp*cosgam)>>T_P)+((cosalp*((sinbet*singam)>>T_P))>>T_P);
  145.  T_sz1=((cosbet*singam)>>T_P);
  146.  T_sx2=-((sinalp*cosbet)>>T_P);
  147.  T_sy2=((cosalp*cosbet)>>T_P);
  148.  T_sz2=-sinbet;
  149.  T_sx3=-((cosalp*singam)>>T_P)-((sinalp*((sinbet*cosgam)>>T_P))>>T_P);
  150.  T_sy3=((cosalp*((sinbet*cosgam)>>T_P))>>T_P)-((sinalp*singam)>>T_P);
  151.  T_sz3=((cosbet*cosgam)>>T_P);              /* calculating the matrix */
  152. #endif
  153. #if defined(_FLOAT_)
  154.  T_sx1=cosalp*cosgam-sinalp*sinbet*singam;
  155.  T_sy1=sinalp*cosgam+cosalp*sinbet*singam;
  156.  T_sz1=cosbet*singam;
  157.  T_sx2=-sinalp*cosbet;
  158.  T_sy2=cosalp*cosbet;
  159.  T_sz2=-sinbet;
  160.  T_sx3=-cosalp*singam-sinalp*sinbet*cosgam;
  161.  T_sy3=cosalp*sinbet*cosgam-sinalp*singam;
  162.  T_sz3=cosbet*cosgam;                       /* calculating the matrix */
  163. #endif
  164. }
  165. /**********************************************************
  166.  * Rotating the coordinates object to world.              *
  167.  *                                                        *
  168.  *                                       |sx1 sx2 sx3|    *
  169.  * T'=T[S]  where:  [x' y' z'] = [x y z]*|sy1 sy2 sy3|    *
  170.  *                                       |sz1 sz2 sz3|    *
  171.  *                                                        *
  172. **********************************************************/
  173. void T_self_rotation(int *from,register int *to,int length)
  174. {
  175.  register int i,xt,yt,zt;
  176.  for(i=0;i<length;i++)
  177.  {
  178.   xt=*from++; yt=*from++; zt=*from++;
  179. #if defined(_FIXED_)
  180.   *to++=(int)(((T_sx1*xt)>>T_P) + ((T_sy1*yt)>>T_P) + ((T_sz1*zt)>>T_P));
  181.   *to++=(int)(((T_sx2*xt)>>T_P) + ((T_sy2*yt)>>T_P) + ((T_sz2*zt)>>T_P));
  182.   *to++=(int)(((T_sx3*xt)>>T_P) + ((T_sy3*yt)>>T_P) + ((T_sz3*zt)>>T_P));
  183. #endif
  184. #if defined(_FLOAT_)
  185.   *to++=(int)(T_sx1*xt + T_sy1*yt + T_sz1*zt);
  186.   *to++=(int)(T_sx2*xt + T_sy2*yt + T_sz2*zt);
  187.   *to++=(int)(T_sx3*xt + T_sy3*yt + T_sz3*zt);
  188. #endif
  189.   to++; from++;                             /* 4th element */
  190.  }
  191. }
  192. /**********************************************************
  193.  * Constructing rotation matrix world to view:            *
  194.  * (gam-bet-alp), yaw, then pitch then roll sequence.     *
  195.  * gam-yaw, bet-pitch alp-roll. NB! angles assumed to be  *
  196.  * in pseudo-degrees in the range [0,256).                *
  197.  *                                                        *
  198.  *          Y^    Z           x'=z*sin(gam)+x*cos(gam)    *
  199.  *           |   /            y'=y                        *
  200.  *           |  / alp         z'=z*cos(gam)-x*sin(gam)    *
  201.  *          /|<---+                                       *
  202.  *         | |/   |           x"=x'                       *
  203.  *  -------|-+-----------> X  y"=y'*cos(bet)-z'*sin(bet)  *
  204.  *     bet | |   __           z"=y'*sin(bet)+z'*cos(bet)  *
  205.  *         V/|   /| gam                                   *
  206.  *         /----+             x"'=y"*sin(alp)+x"*cos(alp) *
  207.  *        /  |                y"'=y"*cos(alp)-x"*sin(alp) *
  208.  *       /   |                z"'=z"                      *
  209.  *           |                                            *
  210.  *                                                        *
  211.  * SETS: T_wx1,...,T_wz3.                                 *
  212.  * -----                                                  *
  213. **********************************************************/
  214. void T_set_world_rotation(unsigned char alp,unsigned char bet,
  215.                           unsigned char gam)
  216. {
  217. #if defined(_FIXED_)
  218.  T_fixed cosalp,sinalp,cosbet,sinbet,cosgam,singam;
  219. #endif
  220. #if defined(_FLOAT_)
  221.  float cosalp,sinalp,cosbet,sinbet,cosgam,singam;
  222. #endif
  223.  cosalp=T_cos[alp];
  224.  sinalp=T_sin[alp];
  225.  cosbet=T_cos[bet];
  226.  sinbet=T_sin[bet];
  227.  cosgam=T_cos[gam];
  228.  singam=T_sin[gam];                         /* initializing */
  229. #if defined(_FIXED_)
  230.  T_wx1=((singam*((sinbet*sinalp)>>T_P))>>T_P) + ((cosgam*cosalp)>>T_P);
  231.  T_wy1=((cosbet*sinalp)>>T_P);
  232.  T_wz1=((singam*cosalp)>>T_P) - ((cosgam*((sinbet*sinalp)>>T_P))>>T_P);
  233.  T_wx2=((singam*((sinbet*cosalp)>>T_P))>>T_P) - ((cosgam*sinalp)>>T_P);
  234.  T_wy2=((cosbet*cosalp)>>T_P);
  235.  T_wz2=-((cosgam*((sinbet*cosalp)>>T_P))>>T_P) - ((singam*sinalp)>>T_P);
  236.  T_wx3=-((singam*cosbet)>>T_P);
  237.  T_wy3=sinbet;
  238.  T_wz3=((cosgam*cosbet)>>T_P);              /* calculating the matrix */
  239. #endif
  240. #if defined(_FLOAT_)
  241.  T_wx1=singam*sinbet*sinalp + cosgam*cosalp;
  242.  T_wy1=cosbet*sinalp;
  243.  T_wz1=singam*cosalp - cosgam*sinbet*sinalp;
  244.  T_wx2=singam*sinbet*cosalp - cosgam*sinalp;
  245.  T_wy2=cosbet*cosalp;
  246.  T_wz2=-cosgam*sinbet*cosalp - singam*sinalp;
  247.  T_wx3=-singam*cosbet;
  248.  T_wy3=sinbet;
  249.  T_wz3=cosgam*cosbet;                       /* calculating the matrix */
  250. #endif
  251. }
  252. /**********************************************************
  253.  * Rotating the coordinates world to view.                *
  254.  *                                                        *
  255.  *                                       |wx1 wx2 wx3|    *
  256.  * T'=T[W]  where:  [x' y' z'] = [x y z]*|wy1 wy2 wy3|    *
  257.  *                                       |wz1 wz2 wz3|    *
  258.  *                                                        *
  259. **********************************************************/
  260. void T_world_rotation(int *from,register int *to,int length)
  261. {
  262.  register int i,xt,yt,zt;
  263.  for(i=0;i<length;i++)
  264.  {
  265.   xt=*from++; yt=*from++; zt=*from++;
  266. #if defined(_FIXED_)
  267.   *to++=(int)(((T_wx1*xt)>>T_P) + ((T_wy1*yt)>>T_P) + ((T_wz1*zt)>>T_P));
  268.   *to++=(int)(((T_wx2*xt)>>T_P) + ((T_wy2*yt)>>T_P) + ((T_wz2*zt)>>T_P));
  269.   *to++=(int)(((T_wx3*xt)>>T_P) + ((T_wy3*yt)>>T_P) + ((T_wz3*zt)>>T_P));
  270. #endif
  271. #if defined(_FLOAT_)
  272.   *to++=(int)(T_wx1*xt+T_wy1*yt+T_wz1*zt);
  273.   *to++=(int)(T_wx2*xt+T_wy2*yt+T_wz2*zt);
  274.   *to++=(int)(T_wx3*xt+T_wy3*yt+T_wz3*zt);
  275. #endif
  276.   to++; from++;                             /* 4th element */
  277.  }
  278. }
  279. /**********************************************************
  280.  * Cancatinating matrices of object to world rotation,    *
  281.  * some translation (may have been object and camera      *
  282.  * translations combined earlier) and world to view       *
  283.  * rotation to obtain a single transformation from object *
  284.  * to view space.                                         *
  285.  *                                                        *
  286.  *                     [K]=[S][T][W]                      *
  287.  *                                                        *
  288.  *                                                        *
  289.  * SETS: T_kx1,...,T_kz3,T_tx,T_ty,T_tz.                  *
  290.  * -----                                                  *
  291. **********************************************************/
  292. void T_cancatinate_self_world(int addx,int addy,int addz)
  293. {
  294. #if defined(_FIXED_)
  295.  T_cx1=(int)(((T_sx1*T_wx1)>>T_P)+((T_sx2*T_wy1)>>T_P)+((T_sx3*T_wz1)>>T_P));
  296.  T_cx2=(int)(((T_sx1*T_wx2)>>T_P)+((T_sx2*T_wy2)>>T_P)+((T_sx3*T_wz2)>>T_P));
  297.  T_cx3=(int)(((T_sx1*T_wx3)>>T_P)+((T_sx2*T_wy3)>>T_P)+((T_sx3*T_wz3)>>T_P));
  298.  T_cy1=(int)(((T_sy1*T_wx1)>>T_P)+((T_sy2*T_wy1)>>T_P)+((T_sy3*T_wz1)>>T_P));
  299.  T_cy2=(int)(((T_sy1*T_wx2)>>T_P)+((T_sy2*T_wy2)>>T_P)+((T_sy3*T_wz2)>>T_P));
  300.  T_cy3=(int)(((T_sy1*T_wx3)>>T_P)+((T_sy2*T_wy3)>>T_P)+((T_sy3*T_wz3)>>T_P));
  301.  T_cz1=(int)(((T_sz1*T_wx1)>>T_P)+((T_sz2*T_wy1)>>T_P)+((T_sz3*T_wz1)>>T_P));
  302.  T_cz2=(int)(((T_sz1*T_wx2)>>T_P)+((T_sz2*T_wy2)>>T_P)+((T_sz3*T_wz2)>>T_P));
  303.  T_cz3=(int)(((T_sz1*T_wx3)>>T_P)+((T_sz2*T_wy3)>>T_P)+((T_sz3*T_wz3)>>T_P));
  304.  T_tx=(int)(((addx*T_wx1)>>T_P)+((addy*T_wy1)>>T_P)+((addz*T_wz1)>>T_P));
  305.  T_ty=(int)(((addx*T_wx2)>>T_P)+((addy*T_wy2)>>T_P)+((addz*T_wz2)>>T_P));
  306.  T_tz=(int)(((addx*T_wx3)>>T_P)+((addy*T_wy3)>>T_P)+((addz*T_wz3)>>T_P));
  307. #endif
  308. #if defined(_FLOAT_)
  309.  T_cx1=T_sx1*T_wx1+T_sx2*T_wy1+T_sx3*T_wz1;
  310.  T_cx2=T_sx1*T_wx2+T_sx2*T_wy2+T_sx3*T_wz2;
  311.  T_cx3=T_sx1*T_wx3+T_sx2*T_wy3+T_sx3*T_wz3;
  312.  T_cy1=T_sy1*T_wx1+T_sy2*T_wy1+T_sy3*T_wz1;
  313.  T_cy2=T_sy1*T_wx2+T_sy2*T_wy2+T_sy3*T_wz2;
  314.  T_cy3=T_sy1*T_wx3+T_sy2*T_wy3+T_sy3*T_wz3;
  315.  T_cz1=T_sz1*T_wx1+T_sz2*T_wy1+T_sz3*T_wz1;
  316.  T_cz2=T_sz1*T_wx2+T_sz2*T_wy2+T_sz3*T_wz2;
  317.  T_cz3=T_sz1*T_wx3+T_sz2*T_wy3+T_sz3*T_wz3;
  318.  T_tx=addx*T_wx1+addy*T_wy1+addz*T_wz1;
  319.  T_ty=addx*T_wx2+addy*T_wy2+addz*T_wz2;
  320.  T_tz=addx*T_wx3+addy*T_wy3+addz*T_wz3;
  321. #endif
  322. }
  323. /**********************************************************
  324.  * Transforming the coordinates object to view.           *
  325.  *                                                        *
  326.  *                                       |cx1 cx2 cx3 0|  *
  327.  * T'=T[K] where: [x' y' z' 1]=[x y z 1]*|cy1 cy2 cy3 0|  *
  328.  *                                       |cz1 cz2 cz3 0|  *
  329.  *                                       |tx  ty  tz  1|  *
  330.  *                                                        *
  331. **********************************************************/
  332. void T_cancatinated_rotation(int *from,register int *to,int length)
  333. {
  334.  register int i,xt,yt,zt;
  335.  for(i=0;i<length;i++)
  336.  {
  337.   xt=*from++; yt=*from++; zt=*from++;
  338. #if defined(_FIXED_)
  339.   *to++=(int)(((T_cx1*xt)>>T_P)+((T_cy1*yt)>>T_P)+((T_cz1*zt)>>T_P)+T_tx);
  340.   *to++=(int)(((T_cx2*xt)>>T_P)+((T_cy2*yt)>>T_P)+((T_cz2*zt)>>T_P)+T_ty);
  341.   *to++=(int)(((T_cx3*xt)>>T_P)+((T_cy3*yt)>>T_P)+((T_cz3*zt)>>T_P)+T_tz);
  342. #endif
  343. #if defined(_FLOAT_)
  344.   *to++=(int)(T_cx1*xt+T_cy1*yt+T_cz1*zt+T_tx);
  345.   *to++=(int)(T_cx2*xt+T_cy2*yt+T_cz2*zt+T_ty);
  346.   *to++=(int)(T_cx3*xt+T_cy3*yt+T_cz3*zt+T_tz);
  347. #endif
  348.   to++; from++;                             /* 4th element */
  349.  }
  350. }
  351. /**********************************************************
  352.  * Transforming to perspective, the coordinates passed    *
  353.  * are supposed to be both volume, and Z-clipped,         *
  354.  * otherwise division by 0 or overflow may occur.         *
  355.  * Also, performs the screen centre translation.          *
  356.  *                                                        *
  357.  *                *                                       *
  358.  *               /|X                                      *
  359.  *              / |                                       *
  360.  *             /  |                                       *
  361.  *            *   |                                       *
  362.  *           /|X' |       X'      X                       *
  363.  *          / |   |    ------- = ---                      *
  364.  *         /  |   |     focus     Z                       *
  365.  *        *---+---+                                       *
  366.  *        0   ^   Z    X'= X*focus/Z                      *
  367.  *            |                                           *
  368.  *          focus                                         *
  369.  *                                                        *
  370. **********************************************************/
  371. void T_perspective(int *from,register int *to,int dimension,
  372.                    int length,int log_focus
  373.                   )
  374. {
  375.  register int i;
  376.  int rem_dim=dimension-3;                   /* other then X Y Z */
  377.  for(i=0;i<length;i++,from+=rem_dim,to+=rem_dim)
  378.  {                                          /* NB! Z is not changed */
  379.   *to++=(int)((((long)from[0])<<log_focus)/from[2])+HW_SCREEN_X_CENTRE;
  380.   *to++=(int)((((long)from[1])<<log_focus)/from[2])+HW_SCREEN_Y_CENTRE;
  381.                                             /* also translates to screen */
  382. #if defined(_Z_BUFFER_)                     /* centre */
  383.   *to++=T_ONE_OVER_Z_CONST/(long)from[2];   /* 1/Z is left in the tuple */
  384. #endif
  385.   HW_copy_int(from+=3,to,rem_dim);          /* remaining values are passed */
  386.  }
  387. }
  388. /**********************************************************
  389.  * Gaussian elimination.                                  *
  390. **********************************************************/
  391. void T_linear_solve(int ia[T_MAX_MATRIX_SIZE][T_MAX_MATRIX_SIZE],
  392.                     int ib[T_MAX_MATRIX_SIZE][T_MAX_MATRIX_SIZE],
  393.                     int ix[T_MAX_MATRIX_SIZE][T_MAX_MATRIX_SIZE],
  394.                     int n,int m,
  395.                     int log_source_multiplyer,
  396.                     int log_result_multiplyer
  397.                    )
  398. {
  399. #if defined(_FIXED_)
  400.  T_fixed a[T_MAX_MATRIX_SIZE][T_MAX_MATRIX_SIZE];
  401.  T_fixed b[T_MAX_MATRIX_SIZE][T_MAX_MATRIX_SIZE];
  402.  T_fixed x[T_MAX_MATRIX_SIZE][T_MAX_MATRIX_SIZE];
  403.  T_fixed max,tmp,pivot,sum;
  404. #endif
  405. #if defined(_FLOAT_)
  406.  float a[T_MAX_MATRIX_SIZE][T_MAX_MATRIX_SIZE];
  407.  float b[T_MAX_MATRIX_SIZE][T_MAX_MATRIX_SIZE];
  408.  float x[T_MAX_MATRIX_SIZE][T_MAX_MATRIX_SIZE];
  409.  float max,tmp,pivot,sum;
  410. #endif
  411.  int i,j,k,num;
  412.  for(i=0;i<n;i++)                           /* into internal representation */
  413.  {
  414. #if defined(_FIXED_)
  415.   for(j=0;j<n;j++) a[i][j]=(T_fixed)(ia[i][j]>>log_source_multiplyer);
  416.   for(j=0;j<m;j++) b[i][j]=(T_fixed)ib[i][j];
  417. #endif
  418. #if defined(_FLOAT_)
  419.   for(j=0;j<n;j++) a[i][j]=(float)(ia[i][j]>>log_source_multiplyer);
  420.   for(j=0;j<m;j++) b[i][j]=(float)ib[i][j];
  421. #endif
  422.  }
  423.  for(max=0,num=0,i=0;i<n-1;i++)
  424.  {
  425.   for(j=i;j<n;j++)                          /* looking for pivot element */
  426.   {
  427.    if(a[j][i]>=0) pivot=a[j][i]; else pivot=-a[j][i];
  428.    if(pivot>max) { max=pivot; num=j; }
  429.   }
  430.   if(max!=0)                                /* for non zero pivots */
  431.   {
  432.    if(num!=i)                               /* swap is requied */
  433.    {
  434.     for(j=0;j<n;j++) { tmp=a[i][j]; a[i][j]=a[num][j]; a[num][j]=tmp; }
  435.     for(j=0;j<m;j++) { tmp=b[i][j]; b[i][j]=b[num][j]; b[num][j]=tmp; }
  436.    }
  437.    for(j=i+1;j<n;j++)                       /* eliminating all coefs below */
  438.    {
  439.     for(k=i+1;k<n;k++) a[j][k]=a[j][k]-a[i][k]*a[j][i]/a[i][i];
  440.     for(k=0;k<m;k++)   b[j][k]=b[j][k]-b[i][k]*a[j][i]/a[i][i];
  441.    }
  442.   }
  443.  }
  444.  for(i=n-1;i>=0;i--)                        /* reversed direction */
  445.  {
  446.   for(k=0;k<m;k++)
  447.   {
  448. #if defined(_FIXED_)
  449.    for(sum=0,j=i+1;j<n;j++) sum=sum+((a[i][j]*x[j][k])>>log_result_multiplyer);
  450.    x[i][k]=(((b[i][k])-sum)<<log_result_multiplyer)/a[i][i];
  451.    ix[i][k]=(int)x[i][k];                   /* external form for the result */
  452. #endif
  453. #if defined(_FLOAT_)
  454.    for(sum=0,j=i+1;j<n;j++) sum=sum+a[i][j]*x[j][k];
  455.    x[i][k]=(b[i][k]-sum)/a[i][i];           /* external form for the result */
  456.    ix[i][k]=(int)(x[i][k]*(0x1<<log_result_multiplyer));
  457. #endif
  458.   }
  459.  }
  460. }
  461. /**********************************************************/