simgear/math/SGMatrix.hxx

   1 // Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Library General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Library General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Library General Public
  14 // License along with this library; if not, write to the
  15 // Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  16 // Boston, MA  02111-1307, USA.
  17 //
  18
  19 #ifndef SGMatrix_H
  20 #define SGMatrix_H
  21
  22 /// Expression templates for poor programmers ... :)
  23 template<typename T>
  24 struct TransNegRef;
  25
  26 template<typename T>
  27 class SGMatrix;
  28
  29 /// 3D Matrix Class
  30 template<typename T>
  31 class SGMatrix {
  32 public:
  33   enum { nCols = 4, nRows = 4, nEnts = 16 };
  34   typedef T value_type;
  35
  36   /// Default constructor. Does not initialize at all.
  37   /// If you need them zero initialized, use SGMatrix::zeros()
  38   SGMatrix(void)
  39   {
  40     /// Initialize with nans in the debug build, that will guarantee to have
  41     /// a fast uninitialized default constructor in the release but shows up
  42     /// uninitialized values in the debug build very fast ...
  43 #ifndef NDEBUG
  44     for (unsigned i = 0; i < nEnts; ++i)
  45       _data.flat[i] = SGLimits<T>::quiet_NaN();
  46 #endif
  47   }
  48   /// Constructor. Initialize by the content of a plain array,
  49   /// make sure it has at least 16 elements
  50   explicit SGMatrix(const T* data)
  51   { for (unsigned i = 0; i < nEnts; ++i) _data.flat[i] = data[i]; }
  52
  53   /// Constructor, build up a SGMatrix from given elements
  54   SGMatrix(T m00, T m01, T m02, T m03,
  55            T m10, T m11, T m12, T m13,
  56            T m20, T m21, T m22, T m23,
  57            T m30, T m31, T m32, T m33)
  58   {
  59     _data.flat[0] = m00; _data.flat[1] = m10;
  60     _data.flat[2] = m20; _data.flat[3] = m30;
  61     _data.flat[4] = m01; _data.flat[5] = m11;
  62     _data.flat[6] = m21; _data.flat[7] = m31;
  63     _data.flat[8] = m02; _data.flat[9] = m12;
  64     _data.flat[10] = m22; _data.flat[11] = m32;
  65     _data.flat[12] = m03; _data.flat[13] = m13;
  66     _data.flat[14] = m23; _data.flat[15] = m33;
  67   }
  68
  69   /// Constructor, build up a SGMatrix from a translation
  70   SGMatrix(const SGVec3<T>& trans)
  71   { set(trans); }
  72
  73   /// Constructor, build up a SGMatrix from a rotation and a translation
  74   SGMatrix(const SGQuat<T>& quat, const SGVec3<T>& trans)
  75   { set(quat, trans); }
  76   /// Constructor, build up a SGMatrix from a rotation and a translation
  77   SGMatrix(const SGQuat<T>& quat)
  78   { set(quat); }
  79
  80   /// Copy constructor for a transposed negated matrix
  81   SGMatrix(const TransNegRef<T>& tm)
  82   { set(tm); }
  83
  84   /// Set from a tranlation
  85   void set(const SGVec3<T>& trans)
  86   {
  87     _data.flat[0] = 1; _data.flat[4] = 0;
  88     _data.flat[8] = 0; _data.flat[12] = -trans(0);
  89     _data.flat[1] = 0; _data.flat[5] = 1;
  90     _data.flat[9] = 0; _data.flat[13] = -trans(1);
  91     _data.flat[2] = 0; _data.flat[6] = 0;
  92     _data.flat[10] = 1; _data.flat[14] = -trans(2);
  93     _data.flat[3] = 0; _data.flat[7] = 0;
  94     _data.flat[11] = 0; _data.flat[15] = 1;
  95   }
  96
  97   /// Set from a scale/rotation and tranlation
  98   void set(const SGQuat<T>& quat, const SGVec3<T>& trans)
  99   {
 100     T w = quat.w(); T x = quat.x(); T y = quat.y(); T z = quat.z();
 101     T xx = x*x; T yy = y*y; T zz = z*z;
 102     T wx = w*x; T wy = w*y; T wz = w*z;
 103     T xy = x*y; T xz = x*z; T yz = y*z;
 104     _data.flat[0] = 1-2*(yy+zz); _data.flat[1] = 2*(xy-wz);
 105     _data.flat[2] = 2*(xz+wy); _data.flat[3] = 0;
 106     _data.flat[4] = 2*(xy+wz); _data.flat[5] = 1-2*(xx+zz);
 107     _data.flat[6] = 2*(yz-wx); _data.flat[7] = 0;
 108     _data.flat[8] = 2*(xz-wy); _data.flat[9] = 2*(yz+wx);
 109     _data.flat[10] = 1-2*(xx+yy); _data.flat[11] = 0;
 110     // Well, this one is ugly here, as that xform method on the current
 111     // object needs the above data to be already set ...
 112     SGVec3<T> t = xformVec(trans);
 113     _data.flat[12] = -t(0); _data.flat[13] = -t(1);
 114     _data.flat[14] = -t(2); _data.flat[15] = 1;
 115   }
 116   /// Set from a scale/rotation and tranlation
 117   void set(const SGQuat<T>& quat)
 118   {
 119     T w = quat.w(); T x = quat.x(); T y = quat.y(); T z = quat.z();
 120     T xx = x*x; T yy = y*y; T zz = z*z;
 121     T wx = w*x; T wy = w*y; T wz = w*z;
 122     T xy = x*y; T xz = x*z; T yz = y*z;
 123     _data.flat[0] = 1-2*(yy+zz); _data.flat[1] = 2*(xy-wz);
 124     _data.flat[2] = 2*(xz+wy); _data.flat[3] = 0;
 125     _data.flat[4] = 2*(xy+wz); _data.flat[5] = 1-2*(xx+zz);
 126     _data.flat[6] = 2*(yz-wx); _data.flat[7] = 0;
 127     _data.flat[8] = 2*(xz-wy); _data.flat[9] = 2*(yz+wx);
 128     _data.flat[10] = 1-2*(xx+yy); _data.flat[11] = 0;
 129     _data.flat[12] = 0; _data.flat[13] = 0;
 130     _data.flat[14] = 0; _data.flat[15] = 1;
 131   }
 132
 133   /// set from a transposed negated matrix
 134   void set(const TransNegRef<T>& tm)
 135   {
 136     const SGMatrix& m = tm.m;
 137     _data.flat[0] = m(0,0);
 138     _data.flat[1] = m(0,1);
 139     _data.flat[2] = m(0,2);
 140     _data.flat[3] = m(3,0);
 141
 142     _data.flat[4] = m(1,0);
 143     _data.flat[5] = m(1,1);
 144     _data.flat[6] = m(1,2);
 145     _data.flat[7] = m(3,1);
 146
 147     _data.flat[8] = m(2,0);
 148     _data.flat[9] = m(2,1);
 149     _data.flat[10] = m(2,2);
 150     _data.flat[11] = m(3,2);
 151
 152     // Well, this one is ugly here, as that xform method on the current
 153     // object needs the above data to be already set ...
 154     SGVec3<T> t = xformVec(SGVec3<T>(m(0,3), m(1,3), m(2,3)));
 155     _data.flat[12] = -t(0);
 156     _data.flat[13] = -t(1);
 157     _data.flat[14] = -t(2);
 158     _data.flat[15] = m(3,3);
 159   }
 160
 161   /// Access by index, the index is unchecked
 162   const T& operator()(unsigned i, unsigned j) const
 163   { return _data.flat[i + 4*j]; }
 164   /// Access by index, the index is unchecked
 165   T& operator()(unsigned i, unsigned j)
 166   { return _data.flat[i + 4*j]; }
 167
 168   /// Access raw data by index, the index is unchecked
 169   const T& operator[](unsigned i) const
 170   { return _data.flat[i]; }
 171   /// Access by index, the index is unchecked
 172   T& operator[](unsigned i)
 173   { return _data.flat[i]; }
 174
 175   /// Get the data pointer
 176   const T* data(void) const
 177   { return _data.flat; }
 178   /// Get the data pointer
 179   T* data(void)
 180   { return _data.flat; }
 181
 182   /// Readonly interface function to ssg's sgMat4/sgdMat4
 183   const T (&sg(void) const)[4][4]
 184   { return _data.carray; }
 185   /// Interface function to ssg's sgMat4/sgdMat4
 186   T (&sg(void))[4][4]
 187   { return _data.carray; }
 188
 189   /// Inplace addition
 190   SGMatrix& operator+=(const SGMatrix& m)
 191   { for (unsigned i = 0; i < nEnts; ++i) _data.flat[i] += m._data.flat[i]; return *this; }
 192   /// Inplace subtraction
 193   SGMatrix& operator-=(const SGMatrix& m)
 194   { for (unsigned i = 0; i < nEnts; ++i) _data.flat[i] -= m._data.flat[i]; return *this; }
 195   /// Inplace scalar multiplication
 196   template<typename S>
 197   SGMatrix& operator*=(S s)
 198   { for (unsigned i = 0; i < nEnts; ++i) _data.flat[i] *= s; return *this; }
 199   /// Inplace scalar multiplication by 1/s
 200   template<typename S>
 201   SGMatrix& operator/=(S s)
 202   { return operator*=(1/T(s)); }
 203   /// Inplace matrix multiplication, post multiply
 204   SGMatrix& operator*=(const SGMatrix<T>& m2);
 205
 206   SGVec3<T> xformPt(const SGVec3<T>& pt) const
 207   {
 208     SGVec3<T> tpt;
 209     tpt(0) = (*this)(0,3);
 210     tpt(1) = (*this)(1,3);
 211     tpt(2) = (*this)(2,3);
 212     for (unsigned i = 0; i < SGMatrix<T>::nCols-1; ++i) {
 213       T tmp = pt(i);
 214       tpt(0) += tmp*(*this)(0,i);
 215       tpt(1) += tmp*(*this)(1,i);
 216       tpt(2) += tmp*(*this)(2,i);
 217     }
 218     return tpt;
 219   }
 220   SGVec3<T> xformVec(const SGVec3<T>& v) const
 221   {
 222     SGVec3<T> tv;
 223     T tmp = v(0);
 224     tv(0) = tmp*(*this)(0,0);
 225     tv(1) = tmp*(*this)(1,0);
 226     tv(2) = tmp*(*this)(2,0);
 227     for (unsigned i = 1; i < SGMatrix<T>::nCols-1; ++i) {
 228       T tmp = v(i);
 229       tv(0) += tmp*(*this)(0,i);
 230       tv(1) += tmp*(*this)(1,i);
 231       tv(2) += tmp*(*this)(2,i);
 232     }
 233     return tv;
 234   }
 235
 236   /// Return an all zero matrix
 237   static SGMatrix zeros(void)
 238   { return SGMatrix(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); }
 239
 240   /// Return a unit matrix
 241   static SGMatrix unit(void)
 242   { return SGMatrix(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1); }
 243
 244 private:
 245   /// Required to make that alias safe.
 246   union Data {
 247     T flat[16];
 248     T carray[4][4];
 249   };
 250
 251   /// The actual data, the matrix is stored in column major order,
 252   /// that matches the storage format of OpenGL
 253   Data _data;
 254 };
 255
 256 /// Class to distinguish between a matrix and the matrix with a transposed
 257 /// rotational part and a negated translational part
 258 template<typename T>
 259 struct TransNegRef {
 260   TransNegRef(const SGMatrix<T>& _m) : m(_m) {}
 261   const SGMatrix<T>& m;
 262 };
 263
 264 /// Unary +, do nothing ...
 265 template<typename T>
 266 inline
 267 const SGMatrix<T>&
 268 operator+(const SGMatrix<T>& m)
 269 { return m; }
 270
 271 /// Unary -, do nearly nothing
 272 template<typename T>
 273 inline
 274 SGMatrix<T>
 275 operator-(const SGMatrix<T>& m)
 276 {
 277   SGMatrix<T> ret;
 278   for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
 279     ret[i] = -m[i];
 280   return ret;
 281 }
 282
 283 /// Binary +
 284 template<typename T>
 285 inline
 286 SGMatrix<T>
 287 operator+(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
 288 {
 289   SGMatrix<T> ret;
 290   for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
 291     ret[i] = m1[i] + m2[i];
 292   return ret;
 293 }
 294
 295 /// Binary -
 296 template<typename T>
 297 inline
 298 SGMatrix<T>
 299 operator-(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
 300 {
 301   SGMatrix<T> ret;
 302   for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
 303     ret[i] = m1[i] - m2[i];
 304   return ret;
 305 }
 306
 307 /// Scalar multiplication
 308 template<typename S, typename T>
 309 inline
 310 SGMatrix<T>
 311 operator*(S s, const SGMatrix<T>& m)
 312 {
 313   SGMatrix<T> ret;
 314   for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
 315     ret[i] = s*m[i];
 316   return ret;
 317 }
 318
 319 /// Scalar multiplication
 320 template<typename S, typename T>
 321 inline
 322 SGMatrix<T>
 323 operator*(const SGMatrix<T>& m, S s)
 324 {
 325   SGMatrix<T> ret;
 326   for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
 327     ret[i] = s*m[i];
 328   return ret;
 329 }
 330
 331 /// Vector multiplication
 332 template<typename T>
 333 inline
 334 SGVec4<T>
 335 operator*(const SGMatrix<T>& m, const SGVec4<T>& v)
 336 {
 337   SGVec4<T> mv;
 338   T tmp = v(0);
 339   mv(0) = tmp*m(0,0);
 340   mv(1) = tmp*m(1,0);
 341   mv(2) = tmp*m(2,0);
 342   mv(3) = tmp*m(3,0);
 343   for (unsigned i = 1; i < SGMatrix<T>::nCols; ++i) {
 344     T tmp = v(i);
 345     mv(0) += tmp*m(0,i);
 346     mv(1) += tmp*m(1,i);
 347     mv(2) += tmp*m(2,i);
 348     mv(3) += tmp*m(3,i);
 349   }
 350   return mv;
 351 }
 352
 353 /// Vector multiplication
 354 template<typename T>
 355 inline
 356 SGVec4<T>
 357 operator*(const TransNegRef<T>& tm, const SGVec4<T>& v)
 358 {
 359   const SGMatrix<T>& m = tm.m;
 360   SGVec4<T> mv;
 361   SGVec3<T> v2;
 362   T tmp = v(3);
 363   mv(0) = v2(0) = -tmp*m(0,3);
 364   mv(1) = v2(1) = -tmp*m(1,3);
 365   mv(2) = v2(2) = -tmp*m(2,3);
 366   mv(3) = tmp*m(3,3);
 367   for (unsigned i = 0; i < SGMatrix<T>::nCols - 1; ++i) {
 368     T tmp = v(i) + v2(i);
 369     mv(0) += tmp*m(i,0);
 370     mv(1) += tmp*m(i,1);
 371     mv(2) += tmp*m(i,2);
 372     mv(3) += tmp*m(3,i);
 373   }
 374   return mv;
 375 }
 376
 377 /// Matrix multiplication
 378 template<typename T>
 379 inline
 380 SGMatrix<T>
 381 operator*(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
 382 {
 383   SGMatrix<T> m;
 384   for (unsigned j = 0; j < SGMatrix<T>::nCols; ++j) {
 385     T tmp = m2(0,j);
 386     m(0,j) = tmp*m1(0,0);
 387     m(1,j) = tmp*m1(1,0);
 388     m(2,j) = tmp*m1(2,0);
 389     m(3,j) = tmp*m1(3,0);
 390     for (unsigned i = 1; i < SGMatrix<T>::nCols; ++i) {
 391       T tmp = m2(i,j);
 392       m(0,j) += tmp*m1(0,i);
 393       m(1,j) += tmp*m1(1,i);
 394       m(2,j) += tmp*m1(2,i);
 395       m(3,j) += tmp*m1(3,i);
 396     }
 397   }
 398   return m;
 399 }
 400
 401 /// Inplace matrix multiplication, post multiply
 402 template<typename T>
 403 inline
 404 SGMatrix<T>&
 405 SGMatrix<T>::operator*=(const SGMatrix<T>& m2)
 406 { (*this) = operator*(*this, m2); return *this; }
 407
 408 /// Return a reference to the matrix typed to make sure we use the transposed
 409 /// negated matrix
 410 template<typename T>
 411 inline
 412 TransNegRef<T>
 413 transNeg(const SGMatrix<T>& m)
 414 { return TransNegRef<T>(m); }
 415
 416 /// Compute the inverse if the matrix src. Store the result in dst.
 417 /// Return if the matrix is nonsingular. If it is singular dst contains
 418 /// undefined values
 419 template<typename T>
 420 inline
 421 bool
 422 invert(SGMatrix<T>& dst, const SGMatrix<T>& src)
 423 {
 424   // Do a LU decomposition with row pivoting and solve into dst
 425   SGMatrix<T> tmp = src;
 426   dst = SGMatrix<T>::unit();
 427
 428   for (unsigned i = 0; i < 4; ++i) {
 429     T val = tmp(i,i);
 430     unsigned ind = i;
 431
 432     // Find the row with the maximum value in the i-th colum
 433     for (unsigned j = i + 1; j < 4; ++j) {
 434       if (fabs(tmp(j, i)) > fabs(val)) {
 435         ind = j;
 436         val = tmp(j, i);
 437       }
 438     }
 439
 440     // Do row pivoting
 441     if (ind != i) {
 442       for (unsigned j = 0; j < 4; ++j) {
 443         T t;
 444         t = dst(i,j); dst(i,j) = dst(ind,j); dst(ind,j) = t;
 445         t = tmp(i,j); tmp(i,j) = tmp(ind,j); tmp(ind,j) = t;
 446       }
 447     }
 448
 449     // Check for singularity
 450     if (fabs(val) <= SGLimits<T>::min())
 451       return false;
 452
 453     T ival = 1/val;
 454     for (unsigned j = 0; j < 4; ++j) {
 455       tmp(i,j) *= ival;
 456       dst(i,j) *= ival;
 457     }
 458
 459     for (unsigned j = 0; j < 4; ++j) {
 460       if (j == i)
 461         continue;
 462
 463       val = tmp(j,i);
 464       for (unsigned k = 0; k < 4; ++k) {
 465         tmp(j,k) -= tmp(i,k) * val;
 466         dst(j,k) -= dst(i,k) * val;
 467       }
 468     }
 469   }
 470   return true;
 471 }
 472
 473 /// The 1-norm of the matrix, this is the largest column sum of
 474 /// the absolute values of A.
 475 template<typename T>
 476 inline
 477 T
 478 norm1(const SGMatrix<T>& m)
 479 {
 480   T nrm = 0;
 481   for (unsigned i = 0; i < SGMatrix<T>::nRows; ++i) {
 482     T sum = fabs(m(i, 0)) + fabs(m(i, 1)) + fabs(m(i, 2)) + fabs(m(i, 3));
 483     if (nrm < sum)
 484       nrm = sum;
 485   }
 486   return nrm;
 487 }
 488
 489 /// The inf-norm of the matrix, this is the largest row sum of
 490 /// the absolute values of A.
 491 template<typename T>
 492 inline
 493 T
 494 normInf(const SGMatrix<T>& m)
 495 {
 496   T nrm = 0;
 497   for (unsigned i = 0; i < SGMatrix<T>::nCols; ++i) {
 498     T sum = fabs(m(0, i)) + fabs(m(1, i)) + fabs(m(2, i)) + fabs(m(3, i));
 499     if (nrm < sum)
 500       nrm = sum;
 501   }
 502   return nrm;
 503 }
 504
 505 /// Return true if exactly the same
 506 template<typename T>
 507 inline
 508 bool
 509 operator==(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
 510 {
 511   for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
 512     if (m1[i] != m2[i])
 513       return false;
 514   return true;
 515 }
 516
 517 /// Return true if not exactly the same
 518 template<typename T>
 519 inline
 520 bool
 521 operator!=(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
 522 { return ! (m1 == m2); }
 523
 524 /// Return true if equal to the relative tolerance tol
 525 template<typename T>
 526 inline
 527 bool
 528 equivalent(const SGMatrix<T>& m1, const SGMatrix<T>& m2, T rtol, T atol)
 529 { return norm1(m1 - m2) < rtol*(norm1(m1) + norm1(m2)) + atol; }
 530
 531 /// Return true if equal to the relative tolerance tol
 532 template<typename T>
 533 inline
 534 bool
 535 equivalent(const SGMatrix<T>& m1, const SGMatrix<T>& m2, T rtol)
 536 { return norm1(m1 - m2) < rtol*(norm1(m1) + norm1(m2)); }
 537
 538 /// Return true if about equal to roundoff of the underlying type
 539 template<typename T>
 540 inline
 541 bool
 542 equivalent(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
 543 {
 544   T tol = 100*SGLimits<T>::epsilon();
 545   return equivalent(m1, m2, tol, tol);
 546 }
 547
 548 #ifndef NDEBUG
 549 template<typename T>
 550 inline
 551 bool
 552 isNaN(const SGMatrix<T>& m)
 553 {
 554   for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i) {
 555     if (SGMisc<T>::isNaN(m[i]))
 556       return true;
 557   }
 558   return false;
 559 }
 560 #endif
 561
 562 /// Output to an ostream
 563 template<typename char_type, typename traits_type, typename T>
 564 inline
 565 std::basic_ostream<char_type, traits_type>&
 566 operator<<(std::basic_ostream<char_type, traits_type>& s, const SGMatrix<T>& m)
 567 {
 568   s << "[ " << m(0,0) << ", " << m(0,1) << ", " << m(0,2) << ", " << m(0,3) << "\n";
 569   s << "  " << m(1,0) << ", " << m(1,1) << ", " << m(1,2) << ", " << m(1,3) << "\n";
 570   s << "  " << m(2,0) << ", " << m(2,1) << ", " << m(2,2) << ", " << m(2,3) << "\n";
 571   s << "  " << m(3,0) << ", " << m(3,1) << ", " << m(3,2) << ", " << m(3,3) << " ]";
 572   return s;
 573 }
 574
 575 /// Two classes doing actually the same on different types
 576 typedef SGMatrix<float> SGMatrixf;
 577 typedef SGMatrix<double> SGMatrixd;
 578
 579 inline
 580 SGMatrixf
 581 toMatrixf(const SGMatrixd& m)
 582 {
 583   return SGMatrixf((float)m(0,0), (float)m(0,1), (float)m(0,2), (float)m(0,3),
 584                    (float)m(1,0), (float)m(1,1), (float)m(1,2), (float)m(1,3),
 585                    (float)m(3,0), (float)m(2,1), (float)m(2,2), (float)m(2,3),
 586                    (float)m(4,0), (float)m(4,1), (float)m(4,2), (float)m(4,3));
 587 }
 588
 589 inline
 590 SGMatrixd
 591 toMatrixd(const SGMatrixf& m)
 592 {
 593   return SGMatrixd(m(0,0), m(0,1), m(0,2), m(0,3),
 594                    m(1,0), m(1,1), m(1,2), m(1,3),
 595                    m(3,0), m(2,1), m(2,2), m(2,3),
 596                    m(4,0), m(4,1), m(4,2), m(4,3));
 597 }
 598
 599 #endif