Open Dynamics Engine: odemath.h Source File

00001 /*************************************************************************
00002  *                                                                       *
00003  * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith.       *
00004  * All rights reserved.  Email: russ@q12.org   Web: www.q12.org          *
00005  *                                                                       *
00006  * This library is free software; you can redistribute it and/or         *
00007  * modify it under the terms of EITHER:                                  *
00008  *   (1) The GNU Lesser General Public License as published by the Free  *
00009  *       Software Foundation; either version 2.1 of the License, or (at  *
00010  *       your option) any later version. The text of the GNU Lesser      *
00011  *       General Public License is included with this library in the     *
00012  *       file LICENSE.TXT.                                               *
00013  *   (2) The BSD-style license that is included with this library in     *
00014  *       the file LICENSE-BSD.TXT.                                       *
00015  *                                                                       *
00016  * This library is distributed in the hope that it will be useful,       *
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files    *
00019  * LICENSE.TXT and LICENSE-BSD.TXT for more details.                     *
00020  *                                                                       *
00021  *************************************************************************/
00022 
00023 #ifndef _ODE_ODEMATH_H_
00024 #define _ODE_ODEMATH_H_
00025 
00026 #include <ode/common.h>
00027 
00028 #ifdef __GNUC__
00029 #define PURE_INLINE extern inline
00030 #else
00031 #define PURE_INLINE inline
00032 #endif
00033 
00034 /*
00035  * macro to access elements i,j in an NxM matrix A, independent of the
00036  * matrix storage convention.
00037  */
00038 #define dACCESS33(A,i,j) ((A)[(i)*4+(j)])
00039 
00040 /*
00041  * Macro to test for valid floating point values
00042  */
00043 #define dVALIDVEC3(v) (!(dIsNan(v[0]) || dIsNan(v[1]) || dIsNan(v[2])))
00044 #define dVALIDVEC4(v) (!(dIsNan(v[0]) || dIsNan(v[1]) || dIsNan(v[2]) || dIsNan(v[3])))
00045 #define dVALIDMAT3(m) (!(dIsNan(m[0]) || dIsNan(m[1]) || dIsNan(m[2]) || dIsNan(m[3]) || dIsNan(m[4]) || dIsNan(m[5]) || dIsNan(m[6]) || dIsNan(m[7]) || dIsNan(m[8]) || dIsNan(m[9]) || dIsNan(m[10]) || dIsNan(m[11])))
00046 #define dVALIDMAT4(m) (!(dIsNan(m[0]) || dIsNan(m[1]) || dIsNan(m[2]) || dIsNan(m[3]) || dIsNan(m[4]) || dIsNan(m[5]) || dIsNan(m[6]) || dIsNan(m[7]) || dIsNan(m[8]) || dIsNan(m[9]) || dIsNan(m[10]) || dIsNan(m[11]) || dIsNan(m[12]) || dIsNan(m[13]) || dIsNan(m[14]) || dIsNan(m[15]) ))
00047 
00048 
00049 
00050 /*
00051  * General purpose vector operations with other vectors or constants.
00052  */
00053 
00054 #define dOP(a,op,b,c) \
00055     (a)[0] = ((b)[0]) op ((c)[0]); \
00056     (a)[1] = ((b)[1]) op ((c)[1]); \
00057     (a)[2] = ((b)[2]) op ((c)[2]);
00058 #define dOPC(a,op,b,c) \
00059     (a)[0] = ((b)[0]) op (c); \
00060     (a)[1] = ((b)[1]) op (c); \
00061     (a)[2] = ((b)[2]) op (c);
00062 #define dOPE(a,op,b) \
00063     (a)[0] op ((b)[0]); \
00064     (a)[1] op ((b)[1]); \
00065     (a)[2] op ((b)[2]);
00066 #define dOPEC(a,op,c) \
00067     (a)[0] op (c); \
00068     (a)[1] op (c); \
00069     (a)[2] op (c);
00070 
00077 #define dOPE2(a,op1,b,op2,c) \
00078     (a)[0] op1 ((b)[0]) op2 ((c)[0]); \
00079     (a)[1] op1 ((b)[1]) op2 ((c)[1]); \
00080     (a)[2] op1 ((b)[2]) op2 ((c)[2]);
00081 
00082 /*
00083  * Length, and squared length helpers. dLENGTH returns the length of a dVector3.
00084  * dLENGTHSQUARED return the squared length of a dVector3.
00085  */
00086 
00087 #define dLENGTHSQUARED(a) (((a)[0])*((a)[0]) + ((a)[1])*((a)[1]) + ((a)[2])*((a)[2]))
00088 
00089 #ifdef __cplusplus
00090 
00091 PURE_INLINE dReal dLENGTH (const dReal *a) { return dSqrt(dLENGTHSQUARED(a)); }
00092 
00093 #else
00094 
00095 #define dLENGTH(a) ( dSqrt( ((a)[0])*((a)[0]) + ((a)[1])*((a)[1]) + ((a)[2])*((a)[2]) ) )
00096 
00097 #endif /* __cplusplus */
00098 
00099 
00100 
00101 
00102 
00103 /*
00104  * 3-way dot product. dDOTpq means that elements of `a' and `b' are spaced
00105  * p and q indexes apart respectively. dDOT() means dDOT11.
00106  * in C++ we could use function templates to get all the versions of these
00107  * functions - but on some compilers this will result in sub-optimal code.
00108  */
00109 
00110 #define dDOTpq(a,b,p,q) ((a)[0]*(b)[0] + (a)[p]*(b)[q] + (a)[2*(p)]*(b)[2*(q)])
00111 
00112 #ifdef __cplusplus
00113 
00114 PURE_INLINE dReal dDOT   (const dReal *a, const dReal *b) { return dDOTpq(a,b,1,1); }
00115 PURE_INLINE dReal dDOT13 (const dReal *a, const dReal *b) { return dDOTpq(a,b,1,3); }
00116 PURE_INLINE dReal dDOT31 (const dReal *a, const dReal *b) { return dDOTpq(a,b,3,1); }
00117 PURE_INLINE dReal dDOT33 (const dReal *a, const dReal *b) { return dDOTpq(a,b,3,3); }
00118 PURE_INLINE dReal dDOT14 (const dReal *a, const dReal *b) { return dDOTpq(a,b,1,4); }
00119 PURE_INLINE dReal dDOT41 (const dReal *a, const dReal *b) { return dDOTpq(a,b,4,1); }
00120 PURE_INLINE dReal dDOT44 (const dReal *a, const dReal *b) { return dDOTpq(a,b,4,4); }
00121 
00122 #else
00123 
00124 #define dDOT(a,b)   dDOTpq(a,b,1,1)
00125 #define dDOT13(a,b) dDOTpq(a,b,1,3)
00126 #define dDOT31(a,b) dDOTpq(a,b,3,1)
00127 #define dDOT33(a,b) dDOTpq(a,b,3,3)
00128 #define dDOT14(a,b) dDOTpq(a,b,1,4)
00129 #define dDOT41(a,b) dDOTpq(a,b,4,1)
00130 #define dDOT44(a,b) dDOTpq(a,b,4,4)
00131 
00132 #endif /* __cplusplus */
00133 
00134 
00135 /*
00136  * cross product, set a = b x c. dCROSSpqr means that elements of `a', `b'
00137  * and `c' are spaced p, q and r indexes apart respectively.
00138  * dCROSS() means dCROSS111. `op' is normally `=', but you can set it to
00139  * +=, -= etc to get other effects.
00140  */
00141 
00142 #define dCROSS(a,op,b,c) \
00143 do { \
00144   (a)[0] op ((b)[1]*(c)[2] - (b)[2]*(c)[1]); \
00145   (a)[1] op ((b)[2]*(c)[0] - (b)[0]*(c)[2]); \
00146   (a)[2] op ((b)[0]*(c)[1] - (b)[1]*(c)[0]); \
00147 } while(0)
00148 #define dCROSSpqr(a,op,b,c,p,q,r) \
00149 do { \
00150   (a)[  0] op ((b)[  q]*(c)[2*r] - (b)[2*q]*(c)[  r]); \
00151   (a)[  p] op ((b)[2*q]*(c)[  0] - (b)[  0]*(c)[2*r]); \
00152   (a)[2*p] op ((b)[  0]*(c)[  r] - (b)[  q]*(c)[  0]); \
00153 } while(0)
00154 #define dCROSS114(a,op,b,c) dCROSSpqr(a,op,b,c,1,1,4)
00155 #define dCROSS141(a,op,b,c) dCROSSpqr(a,op,b,c,1,4,1)
00156 #define dCROSS144(a,op,b,c) dCROSSpqr(a,op,b,c,1,4,4)
00157 #define dCROSS411(a,op,b,c) dCROSSpqr(a,op,b,c,4,1,1)
00158 #define dCROSS414(a,op,b,c) dCROSSpqr(a,op,b,c,4,1,4)
00159 #define dCROSS441(a,op,b,c) dCROSSpqr(a,op,b,c,4,4,1)
00160 #define dCROSS444(a,op,b,c) dCROSSpqr(a,op,b,c,4,4,4)
00161 
00162 
00163 /*
00164  * set a 3x3 submatrix of A to a matrix such that submatrix(A)*b = a x b.
00165  * A is stored by rows, and has `skip' elements per row. the matrix is
00166  * assumed to be already zero, so this does not write zero elements!
00167  * if (plus,minus) is (+,-) then a positive version will be written.
00168  * if (plus,minus) is (-,+) then a negative version will be written.
00169  */
00170 
00171 #define dCROSSMAT(A,a,skip,plus,minus) \
00172 do { \
00173   (A)[1] = minus (a)[2]; \
00174   (A)[2] = plus (a)[1]; \
00175   (A)[(skip)+0] = plus (a)[2]; \
00176   (A)[(skip)+2] = minus (a)[0]; \
00177   (A)[2*(skip)+0] = minus (a)[1]; \
00178   (A)[2*(skip)+1] = plus (a)[0]; \
00179 } while(0)
00180 
00181 
00182 /*
00183  * compute the distance between two 3D-vectors
00184  */
00185 
00186 #ifdef __cplusplus
00187 PURE_INLINE dReal dDISTANCE (const dVector3 a, const dVector3 b)
00188    { return dSqrt( (a[0]-b[0])*(a[0]-b[0]) + (a[1]-b[1])*(a[1]-b[1]) + (a[2]-b[2])*(a[2]-b[2]) ); }
00189 #else
00190 #define dDISTANCE(a,b) \
00191    (dSqrt( ((a)[0]-(b)[0])*((a)[0]-(b)[0]) + ((a)[1]-(b)[1])*((a)[1]-(b)[1]) + ((a)[2]-(b)[2])*((a)[2]-(b)[2]) ))
00192 #endif
00193 
00194 
00195 /*
00196  * special case matrix multipication, with operator selection
00197  */
00198 
00199 #define dMULTIPLYOP0_331(A,op,B,C) \
00200 do { \
00201   (A)[0] op dDOT((B),(C)); \
00202   (A)[1] op dDOT((B+4),(C)); \
00203   (A)[2] op dDOT((B+8),(C)); \
00204 } while(0)
00205 #define dMULTIPLYOP1_331(A,op,B,C) \
00206 do { \
00207   (A)[0] op dDOT41((B),(C)); \
00208   (A)[1] op dDOT41((B+1),(C)); \
00209   (A)[2] op dDOT41((B+2),(C)); \
00210 } while(0)
00211 #define dMULTIPLYOP0_133(A,op,B,C) \
00212 do { \
00213   (A)[0] op dDOT14((B),(C)); \
00214   (A)[1] op dDOT14((B),(C+1)); \
00215   (A)[2] op dDOT14((B),(C+2)); \
00216 } while(0)
00217 #define dMULTIPLYOP0_333(A,op,B,C) \
00218 do { \
00219   (A)[0] op dDOT14((B),(C)); \
00220   (A)[1] op dDOT14((B),(C+1)); \
00221   (A)[2] op dDOT14((B),(C+2)); \
00222   (A)[4] op dDOT14((B+4),(C)); \
00223   (A)[5] op dDOT14((B+4),(C+1)); \
00224   (A)[6] op dDOT14((B+4),(C+2)); \
00225   (A)[8] op dDOT14((B+8),(C)); \
00226   (A)[9] op dDOT14((B+8),(C+1)); \
00227   (A)[10] op dDOT14((B+8),(C+2)); \
00228 } while(0)
00229 #define dMULTIPLYOP1_333(A,op,B,C) \
00230 do { \
00231   (A)[0] op dDOT44((B),(C)); \
00232   (A)[1] op dDOT44((B),(C+1)); \
00233   (A)[2] op dDOT44((B),(C+2)); \
00234   (A)[4] op dDOT44((B+1),(C)); \
00235   (A)[5] op dDOT44((B+1),(C+1)); \
00236   (A)[6] op dDOT44((B+1),(C+2)); \
00237   (A)[8] op dDOT44((B+2),(C)); \
00238   (A)[9] op dDOT44((B+2),(C+1)); \
00239   (A)[10] op dDOT44((B+2),(C+2)); \
00240 } while(0)
00241 #define dMULTIPLYOP2_333(A,op,B,C) \
00242 do { \
00243   (A)[0] op dDOT((B),(C)); \
00244   (A)[1] op dDOT((B),(C+4)); \
00245   (A)[2] op dDOT((B),(C+8)); \
00246   (A)[4] op dDOT((B+4),(C)); \
00247   (A)[5] op dDOT((B+4),(C+4)); \
00248   (A)[6] op dDOT((B+4),(C+8)); \
00249   (A)[8] op dDOT((B+8),(C)); \
00250   (A)[9] op dDOT((B+8),(C+4)); \
00251   (A)[10] op dDOT((B+8),(C+8)); \
00252 } while(0)
00253 
00254 #ifdef __cplusplus
00255 
00256 #define DECL template <class TA, class TB, class TC> PURE_INLINE void
00257 
00258 DECL dMULTIPLY0_331(TA *A, const TB *B, const TC *C) { dMULTIPLYOP0_331(A,=,B,C); }
00259 DECL dMULTIPLY1_331(TA *A, const TB *B, const TC *C) { dMULTIPLYOP1_331(A,=,B,C); }
00260 DECL dMULTIPLY0_133(TA *A, const TB *B, const TC *C) { dMULTIPLYOP0_133(A,=,B,C); }
00261 DECL dMULTIPLY0_333(TA *A, const TB *B, const TC *C) { dMULTIPLYOP0_333(A,=,B,C); }
00262 DECL dMULTIPLY1_333(TA *A, const TB *B, const TC *C) { dMULTIPLYOP1_333(A,=,B,C); }
00263 DECL dMULTIPLY2_333(TA *A, const TB *B, const TC *C) { dMULTIPLYOP2_333(A,=,B,C); }
00264 
00265 DECL dMULTIPLYADD0_331(TA *A, const TB *B, const TC *C) { dMULTIPLYOP0_331(A,+=,B,C); }
00266 DECL dMULTIPLYADD1_331(TA *A, const TB *B, const TC *C) { dMULTIPLYOP1_331(A,+=,B,C); }
00267 DECL dMULTIPLYADD0_133(TA *A, const TB *B, const TC *C) { dMULTIPLYOP0_133(A,+=,B,C); }
00268 DECL dMULTIPLYADD0_333(TA *A, const TB *B, const TC *C) { dMULTIPLYOP0_333(A,+=,B,C); }
00269 DECL dMULTIPLYADD1_333(TA *A, const TB *B, const TC *C) { dMULTIPLYOP1_333(A,+=,B,C); }
00270 DECL dMULTIPLYADD2_333(TA *A, const TB *B, const TC *C) { dMULTIPLYOP2_333(A,+=,B,C); }
00271 
00272 #undef DECL
00273 
00274 #else
00275 
00276 #define dMULTIPLY0_331(A,B,C) dMULTIPLYOP0_331(A,=,B,C)
00277 #define dMULTIPLY1_331(A,B,C) dMULTIPLYOP1_331(A,=,B,C)
00278 #define dMULTIPLY0_133(A,B,C) dMULTIPLYOP0_133(A,=,B,C)
00279 #define dMULTIPLY0_333(A,B,C) dMULTIPLYOP0_333(A,=,B,C)
00280 #define dMULTIPLY1_333(A,B,C) dMULTIPLYOP1_333(A,=,B,C)
00281 #define dMULTIPLY2_333(A,B,C) dMULTIPLYOP2_333(A,=,B,C)
00282 
00283 #define dMULTIPLYADD0_331(A,B,C) dMULTIPLYOP0_331(A,+=,B,C)
00284 #define dMULTIPLYADD1_331(A,B,C) dMULTIPLYOP1_331(A,+=,B,C)
00285 #define dMULTIPLYADD0_133(A,B,C) dMULTIPLYOP0_133(A,+=,B,C)
00286 #define dMULTIPLYADD0_333(A,B,C) dMULTIPLYOP0_333(A,+=,B,C)
00287 #define dMULTIPLYADD1_333(A,B,C) dMULTIPLYOP1_333(A,+=,B,C)
00288 #define dMULTIPLYADD2_333(A,B,C) dMULTIPLYOP2_333(A,+=,B,C)
00289 
00290 #endif
00291 
00292 
00293 #ifdef __cplusplus
00294 extern "C" {
00295 #endif
00296 
00297 /*
00298  * normalize 3x1 and 4x1 vectors (i.e. scale them to unit length)
00299  */
00300 ODE_API int  dSafeNormalize3 (dVector3 a);
00301 ODE_API int  dSafeNormalize4 (dVector4 a);
00302 
00303 // For some reason demo_chain1.c does not understand "inline" keyword.
00304 static __inline void _dNormalize3(dVector3 a)
00305 {
00306    int bNormalizationResult = dSafeNormalize3(a);
00307    dIASSERT(bNormalizationResult);
00308    dVARIABLEUSED(bNormalizationResult);
00309 }
00310 
00311 static __inline void _dNormalize4(dVector4 a)
00312 {
00313    int bNormalizationResult = dSafeNormalize4(a);
00314    dIASSERT(bNormalizationResult);
00315    dVARIABLEUSED(bNormalizationResult);
00316 }
00317 
00318 // For DLL export
00319 ODE_API void dNormalize3 (dVector3 a); // Potentially asserts on zero vec
00320 ODE_API void dNormalize4 (dVector4 a); // Potentially asserts on zero vec
00321 
00322 // For internal use
00323 #define dNormalize3(a) _dNormalize3(a)
00324 #define dNormalize4(a) _dNormalize4(a)
00325 
00326 /*
00327  * given a unit length "normal" vector n, generate vectors p and q vectors
00328  * that are an orthonormal basis for the plane space perpendicular to n.
00329  * i.e. this makes p,q such that n,p,q are all perpendicular to each other.
00330  * q will equal n x p. if n is not unit length then p will be unit length but
00331  * q wont be.
00332  */
00333 
00334 ODE_API void dPlaneSpace (const dVector3 n, dVector3 p, dVector3 q);
00335 
00336 #ifdef __cplusplus
00337 }
00338 #endif
00339 
00340 #endif