ab2a30b49f231f192db97ac37eb9ce18202c9a63
[blender.git] / intern / elbeem / intern / ntl_matrices.h
1 /** \file elbeem/intern/ntl_matrices.h
2  *  \ingroup elbeem
3  */
4
5 /******************************************************************************
6  *
7  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
8  * Copyright 2003-2006 Nils Thuerey
9  *
10  * Basic matrix utility include file
11  *
12  *****************************************************************************/
13 #ifndef NTL_MATRICES_H
14
15 #include "ntl_vector3dim.h"
16
17
18 // The basic vector class
19 template<class Scalar>
20 class ntlMatrix4x4
21 {
22 public:
23   // Constructor
24   inline ntlMatrix4x4(void );
25   // Copy-Constructor
26   inline ntlMatrix4x4(const ntlMatrix4x4<Scalar> &v );
27   // construct a vector from one Scalar
28   inline ntlMatrix4x4(Scalar);
29   // construct a vector from three Scalars
30   inline ntlMatrix4x4(Scalar, Scalar, Scalar);
31
32   // Assignment operator
33   inline const ntlMatrix4x4<Scalar>& operator=  (const ntlMatrix4x4<Scalar>& v);
34   // Assignment operator
35   inline const ntlMatrix4x4<Scalar>& operator=  (Scalar s);
36   // Assign and add operator
37   inline const ntlMatrix4x4<Scalar>& operator+= (const ntlMatrix4x4<Scalar>& v);
38   // Assign and add operator
39   inline const ntlMatrix4x4<Scalar>& operator+= (Scalar s);
40   // Assign and sub operator
41   inline const ntlMatrix4x4<Scalar>& operator-= (const ntlMatrix4x4<Scalar>& v);
42   // Assign and sub operator
43   inline const ntlMatrix4x4<Scalar>& operator-= (Scalar s);
44   // Assign and mult operator
45   inline const ntlMatrix4x4<Scalar>& operator*= (const ntlMatrix4x4<Scalar>& v);
46   // Assign and mult operator
47   inline const ntlMatrix4x4<Scalar>& operator*= (Scalar s);
48   // Assign and div operator
49   inline const ntlMatrix4x4<Scalar>& operator/= (const ntlMatrix4x4<Scalar>& v);
50   // Assign and div operator
51   inline const ntlMatrix4x4<Scalar>& operator/= (Scalar s);
52
53   
54   // unary operator
55   inline ntlMatrix4x4<Scalar> operator- () const;
56
57   // binary operator add
58   inline ntlMatrix4x4<Scalar> operator+ (const ntlMatrix4x4<Scalar>&) const;
59   // binary operator add
60   inline ntlMatrix4x4<Scalar> operator+ (Scalar) const;
61   // binary operator sub
62   inline ntlMatrix4x4<Scalar> operator- (const ntlMatrix4x4<Scalar>&) const;
63   // binary operator sub
64   inline ntlMatrix4x4<Scalar> operator- (Scalar) const;
65   // binary operator mult
66   inline ntlMatrix4x4<Scalar> operator* (const ntlMatrix4x4<Scalar>&) const;
67   // binary operator mult
68   inline ntlVector3Dim<Scalar> operator* (const ntlVector3Dim<Scalar>&) const;
69   // binary operator mult
70   inline ntlMatrix4x4<Scalar> operator* (Scalar) const;
71   // binary operator div
72   inline ntlMatrix4x4<Scalar> operator/ (Scalar) const;
73
74         // init function
75         //! init identity matrix
76         inline void initId();
77         //! init rotation matrix
78         inline void initTranslation(Scalar x, Scalar y, Scalar z);
79         //! init rotation matrix
80         inline void initRotationX(Scalar rot);
81         inline void initRotationY(Scalar rot);
82         inline void initRotationZ(Scalar rot);
83         inline void initRotationXYZ(Scalar rotx,Scalar roty, Scalar rotz);
84         //! init scaling matrix
85         inline void initScaling(Scalar scale);
86         inline void initScaling(Scalar x, Scalar y, Scalar z);
87         //! from 16 value array (init id if all 0)
88         inline void initArrayCheck(Scalar *array);
89
90         //! decompose matrix
91         void decompose(ntlVector3Dim<Scalar> &trans,ntlVector3Dim<Scalar> &scale,ntlVector3Dim<Scalar> &rot,ntlVector3Dim<Scalar> &shear);
92
93         //! public to avoid [][] operators
94   Scalar value[4][4];  //< Storage of vector values
95
96
97 protected:
98
99 };
100
101
102
103 //------------------------------------------------------------------------------
104 // TYPEDEFS
105 //------------------------------------------------------------------------------
106
107 // a 3D vector for graphics output, typically float?
108 //typedef ntlMatrix4x4<float>  ntlVec3Gfx; 
109
110 //typedef ntlMatrix4x4<double>  ntlMat4d; 
111 typedef ntlMatrix4x4<double>  ntlMat4d; 
112
113 // a 3D vector with single precision
114 typedef ntlMatrix4x4<float>   ntlMat4f; 
115
116 // a 3D vector with grafix precision
117 typedef ntlMatrix4x4<gfxReal>   ntlMat4Gfx;
118
119 // a 3D integer vector
120 typedef ntlMatrix4x4<int>     ntlMat4i; 
121
122
123
124
125
126 //------------------------------------------------------------------------------
127 // STREAM FUNCTIONS
128 //------------------------------------------------------------------------------
129
130
131
132 /*************************************************************************
133   Outputs the object in human readable form using the format
134   [x,y,z]
135   */
136 template<class Scalar>
137 std::ostream&
138 operator<<( std::ostream& os, const ntlMatrix4x4<Scalar>& m )
139 {
140         for(int i=0; i<4; i++) {
141         os << '|' << m.value[i][0] << ", " << m.value[i][1] << ", " << m.value[i][2] << ", " << m.value[i][3] << '|';
142         }
143   return os;
144 }
145
146
147
148 /*************************************************************************
149   Reads the contents of the object from a stream using the same format
150   as the output operator.
151   */
152 template<class Scalar>
153 std::istream&
154 operator>>( std::istream& is, ntlMatrix4x4<Scalar>& m )
155 {
156   char c;
157   char dummy[3];
158   
159         for(int i=0; i<4; i++) {
160         is >> c >> m.value[i][0] >> dummy >> m.value[i][1] >> dummy >> m.value[i][2] >> dummy >> m.value[i][3] >> c;
161         }
162   return is;
163 }
164
165
166 //------------------------------------------------------------------------------
167 // VECTOR inline FUNCTIONS
168 //------------------------------------------------------------------------------
169
170
171
172 /*************************************************************************
173   Constructor.
174   */
175 template<class Scalar>
176 inline ntlMatrix4x4<Scalar>::ntlMatrix4x4( void )
177 {
178 #ifdef MATRIX_INIT_ZERO
179         for(int i=0; i<4; i++) {
180                 for(int j=0; j<4; j++) {
181                 value[i][j] = 0.0;
182                 }
183         }
184 #endif
185 }
186
187
188
189 /*************************************************************************
190   Copy-Constructor.
191   */
192 template<class Scalar>
193 inline ntlMatrix4x4<Scalar>::ntlMatrix4x4( const ntlMatrix4x4<Scalar> &v )
194 {
195   value[0][0] = v.value[0][0]; value[0][1] = v.value[0][1]; value[0][2] = v.value[0][2]; value[0][3] = v.value[0][3];
196   value[1][0] = v.value[1][0]; value[1][1] = v.value[1][1]; value[1][2] = v.value[1][2]; value[1][3] = v.value[1][3];
197   value[2][0] = v.value[2][0]; value[2][1] = v.value[2][1]; value[2][2] = v.value[2][2]; value[2][3] = v.value[2][3];
198   value[3][0] = v.value[3][0]; value[3][1] = v.value[3][1]; value[3][2] = v.value[3][2]; value[3][3] = v.value[3][3];
199 }
200
201
202
203 /*************************************************************************
204   Constructor for a vector from a single Scalar. All components of
205   the vector get the same value.
206   \param s The value to set
207   \return The new vector
208   */
209 template<class Scalar>
210 inline ntlMatrix4x4<Scalar>::ntlMatrix4x4(Scalar s )
211 {
212         for(int i=0; i<4; i++) {
213                 for(int j=0; j<4; j++) {
214                 value[i][j] = s;
215                 }
216         }
217 }
218
219
220
221 /*************************************************************************
222   Copy a ntlMatrix4x4 componentwise.
223   \param v vector with values to be copied
224   \return Reference to self
225   */
226 template<class Scalar>
227 inline const ntlMatrix4x4<Scalar>&
228 ntlMatrix4x4<Scalar>::operator=( const ntlMatrix4x4<Scalar> &v )
229 {
230   value[0][0] = v.value[0][0]; value[0][1] = v.value[0][1]; value[0][2] = v.value[0][2]; value[0][3] = v.value[0][3];
231   value[1][0] = v.value[1][0]; value[1][1] = v.value[1][1]; value[1][2] = v.value[1][2]; value[1][3] = v.value[1][3];
232   value[2][0] = v.value[2][0]; value[2][1] = v.value[2][1]; value[2][2] = v.value[2][2]; value[2][3] = v.value[2][3];
233   value[3][0] = v.value[3][0]; value[3][1] = v.value[3][1]; value[3][2] = v.value[3][2]; value[3][3] = v.value[3][3];
234   return *this;
235 }
236
237
238 /*************************************************************************
239   Copy a Scalar to each component.
240   \param s The value to copy
241   \return Reference to self
242   */
243 template<class Scalar>
244 inline const ntlMatrix4x4<Scalar>&
245 ntlMatrix4x4<Scalar>::operator=(Scalar s)
246 {
247         for(int i=0; i<4; i++) {
248                 for(int j=0; j<4; j++) {
249                 value[i][j] = s;
250                 }
251         }
252   return *this;
253 }
254
255
256 /*************************************************************************
257   Add another ntlMatrix4x4 componentwise.
258   \param v vector with values to be added
259   \return Reference to self
260   */
261 template<class Scalar>
262 inline const ntlMatrix4x4<Scalar>&
263 ntlMatrix4x4<Scalar>::operator+=( const ntlMatrix4x4<Scalar> &v )
264 {
265   value[0][0] += v.value[0][0]; value[0][1] += v.value[0][1]; value[0][2] += v.value[0][2]; value[0][3] += v.value[0][3];
266   value[1][0] += v.value[1][0]; value[1][1] += v.value[1][1]; value[1][2] += v.value[1][2]; value[1][3] += v.value[1][3];
267   value[2][0] += v.value[2][0]; value[2][1] += v.value[2][1]; value[2][2] += v.value[2][2]; value[2][3] += v.value[2][3];
268   value[3][0] += v.value[3][0]; value[3][1] += v.value[3][1]; value[3][2] += v.value[3][2]; value[3][3] += v.value[3][3];
269   return *this;
270 }
271
272
273 /*************************************************************************
274   Add a Scalar value to each component.
275   \param s Value to add
276   \return Reference to self
277   */
278 template<class Scalar>
279 inline const ntlMatrix4x4<Scalar>&
280 ntlMatrix4x4<Scalar>::operator+=(Scalar s)
281 {
282         for(int i=0; i<4; i++) {
283                 for(int j=0; j<4; j++) {
284                 value[i][j] += s;
285                 }
286         }
287   return *this;
288 }
289
290
291 /*************************************************************************
292   Subtract another vector componentwise.
293   \param v vector of values to subtract
294   \return Reference to self
295   */
296 template<class Scalar>
297 inline const ntlMatrix4x4<Scalar>&
298 ntlMatrix4x4<Scalar>::operator-=( const ntlMatrix4x4<Scalar> &v )
299 {
300   value[0][0] -= v.value[0][0]; value[0][1] -= v.value[0][1]; value[0][2] -= v.value[0][2]; value[0][3] -= v.value[0][3];
301   value[1][0] -= v.value[1][0]; value[1][1] -= v.value[1][1]; value[1][2] -= v.value[1][2]; value[1][3] -= v.value[1][3];
302   value[2][0] -= v.value[2][0]; value[2][1] -= v.value[2][1]; value[2][2] -= v.value[2][2]; value[2][3] -= v.value[2][3];
303   value[3][0] -= v.value[3][0]; value[3][1] -= v.value[3][1]; value[3][2] -= v.value[3][2]; value[3][3] -= v.value[3][3];
304   return *this;
305 }
306
307
308 /*************************************************************************
309   Subtract a Scalar value from each component.
310   \param s Value to subtract
311   \return Reference to self
312   */
313 template<class Scalar>
314 inline const ntlMatrix4x4<Scalar>&
315 ntlMatrix4x4<Scalar>::operator-=(Scalar s)
316 {
317         for(int i=0; i<4; i++) {
318                 for(int j=0; j<4; j++) {
319                 value[i][j] -= s;
320                 }
321         }
322   return *this;
323 }
324
325
326 /*************************************************************************
327   Multiply with another vector componentwise.
328   \param v vector of values to multiply with
329   \return Reference to self
330   */
331 template<class Scalar>
332 inline const ntlMatrix4x4<Scalar>&
333 ntlMatrix4x4<Scalar>::operator*=( const ntlMatrix4x4<Scalar> &v )
334 {
335         ntlMatrix4x4<Scalar> nv(0.0);
336         for(int i=0; i<4; i++) {
337                 for(int j=0; j<4; j++) {
338
339                         for(int k=0;k<4;k++)
340                                 nv.value[i][j] += (value[i][k] * v.value[k][j]);
341                 }
342         }
343   *this = nv;
344   return *this;
345 }
346
347
348 /*************************************************************************
349   Multiply each component with a Scalar value.
350   \param s Value to multiply with
351   \return Reference to self
352   */
353 template<class Scalar>
354 inline const ntlMatrix4x4<Scalar>&
355 ntlMatrix4x4<Scalar>::operator*=(Scalar s)
356 {
357         for(int i=0; i<4; i++) {
358                 for(int j=0; j<4; j++) {
359                 value[i][j] *= s;
360                 }
361         }
362   return *this;
363 }
364
365
366
367 /*************************************************************************
368   Divide each component by a Scalar value.
369   \param s Value to divide by
370   \return Reference to self
371   */
372 template<class Scalar>
373 inline const ntlMatrix4x4<Scalar>&
374 ntlMatrix4x4<Scalar>::operator/=(Scalar s)
375 {
376         for(int i=0; i<4; i++) {
377                 for(int j=0; j<4; j++) {
378                 value[i][j] /= s;
379                 }
380         }
381   return *this;
382 }
383
384
385 //------------------------------------------------------------------------------
386 // unary operators
387 //------------------------------------------------------------------------------
388
389
390 /*************************************************************************
391   Build componentwise the negative this vector.
392   \return The new (negative) vector
393   */
394 template<class Scalar>
395 inline ntlMatrix4x4<Scalar>
396 ntlMatrix4x4<Scalar>::operator-() const
397 {
398         ntlMatrix4x4<Scalar> nv;
399         for(int i=0; i<4; i++) {
400                 for(int j=0; j<4; j++) {
401                 nv[i][j] = -value[i][j];
402                 }
403         }
404   return nv;
405 }
406
407
408
409 //------------------------------------------------------------------------------
410 // binary operators
411 //------------------------------------------------------------------------------
412
413
414 /*************************************************************************
415   Build a vector with another vector added componentwise.
416   \param v The second vector to add
417   \return The sum vector
418   */
419 template<class Scalar>
420 inline ntlMatrix4x4<Scalar>
421 ntlMatrix4x4<Scalar>::operator+( const ntlMatrix4x4<Scalar> &v ) const
422 {
423         ntlMatrix4x4<Scalar> nv;
424         for(int i=0; i<4; i++) {
425                 for(int j=0; j<4; j++) {
426                 nv[i][j] = value[i][j] + v.value[i][j];
427                 }
428         }
429   return nv;
430 }
431
432
433 /*************************************************************************
434   Build a vector with a Scalar value added to each component.
435   \param s The Scalar value to add
436   \return The sum vector
437   */
438 template<class Scalar>
439 inline ntlMatrix4x4<Scalar>
440 ntlMatrix4x4<Scalar>::operator+(Scalar s) const
441 {
442         ntlMatrix4x4<Scalar> nv;
443         for(int i=0; i<4; i++) {
444                 for(int j=0; j<4; j++) {
445                 nv[i][j] = value[i][j] + s;
446                 }
447         }
448   return nv;
449 }
450
451
452 /*************************************************************************
453   Build a vector with another vector subtracted componentwise.
454   \param v The second vector to subtract
455   \return The difference vector
456   */
457 template<class Scalar>
458 inline ntlMatrix4x4<Scalar>
459 ntlMatrix4x4<Scalar>::operator-( const ntlMatrix4x4<Scalar> &v ) const
460 {
461         ntlMatrix4x4<Scalar> nv;
462         for(int i=0; i<4; i++) {
463                 for(int j=0; j<4; j++) {
464                 nv[i][j] = value[i][j] - v.value[i][j];
465                 }
466         }
467   return nv;
468 }
469
470
471 /*************************************************************************
472   Build a vector with a Scalar value subtracted componentwise.
473   \param s The Scalar value to subtract
474   \return The difference vector
475   */
476 template<class Scalar>
477 inline ntlMatrix4x4<Scalar>
478 ntlMatrix4x4<Scalar>::operator-(Scalar s ) const
479 {
480         ntlMatrix4x4<Scalar> nv;
481         for(int i=0; i<4; i++) {
482                 for(int j=0; j<4; j++) {
483                 nv[i][j] = value[i][j] - s;
484                 }
485         }
486   return nv;
487 }
488
489
490
491 /*************************************************************************
492   Build a ntlMatrix4x4 with a Scalar value multiplied to each component.
493   \param s The Scalar value to multiply with
494   \return The product vector
495   */
496 template<class Scalar>
497 inline ntlMatrix4x4<Scalar>
498 ntlMatrix4x4<Scalar>::operator*(Scalar s) const
499 {
500         ntlMatrix4x4<Scalar> nv;
501         for(int i=0; i<4; i++) {
502                 for(int j=0; j<4; j++) {
503                 nv[i][j] = value[i][j] * s;
504                 }
505         }
506   return nv;
507 }
508
509
510
511
512 /*************************************************************************
513   Build a vector divided componentwise by a Scalar value.
514   \param s The Scalar value to divide by
515   \return The ratio vector
516   */
517 template<class Scalar>
518 inline ntlMatrix4x4<Scalar>
519 ntlMatrix4x4<Scalar>::operator/(Scalar s) const
520 {
521         ntlMatrix4x4<Scalar> nv;
522         for(int i=0; i<4; i++) {
523                 for(int j=0; j<4; j++) {
524                 nv[i][j] = value[i][j] / s;
525                 }
526         }
527   return nv;
528 }
529
530
531
532
533
534 /*************************************************************************
535   Build a vector with another vector multiplied by componentwise.
536   \param v The second vector to muliply with
537   \return The product vector
538   */
539 template<class Scalar>
540 inline ntlMatrix4x4<Scalar>
541 ntlMatrix4x4<Scalar>::operator*( const ntlMatrix4x4<Scalar>& v) const
542 {
543         ntlMatrix4x4<Scalar> nv(0.0);
544         for(int i=0; i<4; i++) {
545                 for(int j=0; j<4; j++) {
546
547                         for(int k=0;k<4;k++)
548                                 nv.value[i][j] += (value[i][k] * v.value[k][j]);
549                 }
550         }
551   return nv;
552 }
553
554
555 template<class Scalar>
556 inline ntlVector3Dim<Scalar>
557 ntlMatrix4x4<Scalar>::operator*( const ntlVector3Dim<Scalar>& v) const
558 {
559         ntlVector3Dim<Scalar> nvec(0.0);
560         for(int i=0; i<3; i++) {
561                 for(int j=0; j<3; j++) {
562                         nvec[i] += (v[j] * value[i][j]);
563                 }
564         }
565         // assume normalized w coord
566         for(int i=0; i<3; i++) {
567                 nvec[i] += (1.0 * value[i][3]);
568         }
569   return nvec;
570 }
571
572
573
574 //------------------------------------------------------------------------------
575 // Other helper functions
576 //------------------------------------------------------------------------------
577
578 //! init identity matrix
579 template<class Scalar>
580 inline void ntlMatrix4x4<Scalar>::initId()
581 {
582         (*this) = (Scalar)(0.0);
583         value[0][0] = 
584         value[1][1] = 
585         value[2][2] = 
586         value[3][3] = (Scalar)(1.0);
587 }
588
589 //! init rotation matrix
590 template<class Scalar>
591 inline void ntlMatrix4x4<Scalar>::initTranslation(Scalar x, Scalar y, Scalar z)
592 {
593         //(*this) = (Scalar)(0.0);
594         this->initId();
595         value[0][3] = x;
596         value[1][3] = y;
597         value[2][3] = z;
598 }
599
600 //! init rotation matrix
601 template<class Scalar>
602 inline void 
603 ntlMatrix4x4<Scalar>::initRotationX(Scalar rot)
604 {
605         double drot = (double)(rot/360.0*2.0*M_PI);
606         //? while(drot < 0.0) drot += (M_PI*2.0);
607
608         this->initId();
609         value[1][1] = (Scalar)  cos(drot);
610         value[1][2] = (Scalar)  sin(drot);
611         value[2][1] = (Scalar)(-sin(drot));
612         value[2][2] = (Scalar)  cos(drot);
613 }
614 template<class Scalar>
615 inline void 
616 ntlMatrix4x4<Scalar>::initRotationY(Scalar rot)
617 {
618         double drot = (double)(rot/360.0*2.0*M_PI);
619         //? while(drot < 0.0) drot += (M_PI*2.0);
620
621         this->initId();
622         value[0][0] = (Scalar)  cos(drot);
623         value[0][2] = (Scalar)(-sin(drot));
624         value[2][0] = (Scalar)  sin(drot);
625         value[2][2] = (Scalar)  cos(drot);
626 }
627 template<class Scalar>
628 inline void 
629 ntlMatrix4x4<Scalar>::initRotationZ(Scalar rot)
630 {
631         double drot = (double)(rot/360.0*2.0*M_PI);
632         //? while(drot < 0.0) drot += (M_PI*2.0);
633
634         this->initId();
635         value[0][0] = (Scalar)  cos(drot);
636         value[0][1] = (Scalar)  sin(drot);
637         value[1][0] = (Scalar)(-sin(drot));
638         value[1][1] = (Scalar)  cos(drot);
639 }
640 template<class Scalar>
641 inline void 
642 ntlMatrix4x4<Scalar>::initRotationXYZ( Scalar rotx, Scalar roty, Scalar rotz)
643 {
644         ntlMatrix4x4<Scalar> val;
645         ntlMatrix4x4<Scalar> rot;
646         this->initId();
647
648         // org
649         /*rot.initRotationX(rotx);
650         (*this) *= rot;
651         rot.initRotationY(roty);
652         (*this) *= rot;
653         rot.initRotationZ(rotz);
654         (*this) *= rot;
655         // org */
656
657         // blender
658         rot.initRotationZ(rotz);
659         (*this) *= rot;
660         rot.initRotationY(roty);
661         (*this) *= rot;
662         rot.initRotationX(rotx);
663         (*this) *= rot;
664         // blender */
665 }
666
667 //! init scaling matrix
668 template<class Scalar>
669 inline void 
670 ntlMatrix4x4<Scalar>::initScaling(Scalar scale)
671 {
672         this->initId();
673         value[0][0] = scale;
674         value[1][1] = scale;
675         value[2][2] = scale;
676 }
677 //! init scaling matrix
678 template<class Scalar>
679 inline void 
680 ntlMatrix4x4<Scalar>::initScaling(Scalar x, Scalar y, Scalar z)
681 {
682         this->initId();
683         value[0][0] = x;
684         value[1][1] = y;
685         value[2][2] = z;
686 }
687
688
689 //! from 16 value array (init id if all 0)
690 template<class Scalar>
691 inline void 
692 ntlMatrix4x4<Scalar>::initArrayCheck(Scalar *array)
693 {
694         bool allZero = true;
695         for(int i=0; i<4; i++) {
696                 for(int j=0; j<4; j++) {
697                         value[i][j] = array[i*4+j];
698                         if(array[i*4+j]!=0.0) allZero=false;
699                 }
700         }
701         if(allZero) this->initId();
702 }
703
704 //! decompose matrix
705 template<class Scalar>
706 void 
707 ntlMatrix4x4<Scalar>::decompose(ntlVector3Dim<Scalar> &trans,ntlVector3Dim<Scalar> &scale,ntlVector3Dim<Scalar> &rot,ntlVector3Dim<Scalar> &shear) {
708         ntlVec3Gfx row[3],temp;
709
710         for(int i = 0; i < 3; i++) {
711                 trans[i] = this->value[3][i];
712         }
713
714         for(int i = 0; i < 3; i++) {
715                 row[i][0] = this->value[i][0];
716                 row[i][1] = this->value[i][1];
717                 row[i][2] = this->value[i][2];
718         }
719
720         scale[0] = norm(row[0]);
721         normalize (row[0]);
722
723         shear[0] = dot(row[0], row[1]);
724         row[1][0] = row[1][0] - shear[0]*row[0][0];
725         row[1][1] = row[1][1] - shear[0]*row[0][1];
726         row[1][2] = row[1][2] - shear[0]*row[0][2];
727
728         scale[1] = norm(row[1]);
729         normalize (row[1]);
730
731         if(scale[1] != 0.0)
732                 shear[0] /= scale[1];
733
734         shear[1] = dot(row[0], row[2]);
735         row[2][0] = row[2][0] - shear[1]*row[0][0];
736         row[2][1] = row[2][1] - shear[1]*row[0][1];
737         row[2][2] = row[2][2] - shear[1]*row[0][2];
738
739         shear[2] = dot(row[1], row[2]);
740         row[2][0] = row[2][0] - shear[2]*row[1][0];
741         row[2][1] = row[2][1] - shear[2]*row[1][1];
742         row[2][2] = row[2][2] - shear[2]*row[1][2];
743
744         scale[2] = norm(row[2]);
745         normalize (row[2]);
746
747         if(scale[2] != 0.0) {
748                 shear[1] /= scale[2];
749                 shear[2] /= scale[2];
750         }
751
752         temp = cross(row[1], row[2]);
753         if(dot(row[0], temp) < 0.0) {
754                 for(int i = 0; i < 3; i++) {
755                         scale[i]  *= -1.0;
756                         row[i][0] *= -1.0;
757                         row[i][1] *= -1.0;
758                         row[i][2] *= -1.0;
759                 }
760         }
761
762         if(row[0][2] < -1.0) row[0][2] = -1.0;
763         if(row[0][2] > +1.0) row[0][2] = +1.0;
764
765         rot[1] = asin(-row[0][2]);
766
767         if(fabs(cos(rot[1])) > VECTOR_EPSILON) {
768                 rot[0] = atan2 (row[1][2], row[2][2]);
769                 rot[2] = atan2 (row[0][1], row[0][0]);
770         }
771         else {
772                 rot[0] = atan2 (row[1][0], row[1][1]);
773                 rot[2] = 0.0;
774         }
775
776         rot[0] = (180.0/M_PI)*rot[0];
777         rot[1] = (180.0/M_PI)*rot[1];
778         rot[2] = (180.0/M_PI)*rot[2];
779
780
781 #define NTL_MATRICES_H
782 #endif
783