~ubuntu-branches/ubuntu/maverick/blender/maverick : revision 16

extern/qdune/core/Mathutil.h

// Various math & other functions

#ifndef _MATHUTIL_H

#define _MATHUTIL_H

#include "qdVector.h"

#include "Bound.h"

#define _USE_MATH_DEFINES

#include <cmath>

#include <cassert>

// for mem...()

#include <cstring>

#include "ri.h"

#include "QDRender.h"

__BEGIN_QDRENDER

//------------------------------------------------------------------------------

// various template functions

// returns absolute value of x

template<typename T>

inline T ABS(T x) { return (x < (T)0) ? -x : x; }

// to be used with points/vectors, mostly to test ccw or cw order

// returns twice the (signed) area of triangle defined by vertices (a, b, c)

template<typename T>

inline float signedArea2X(const T &a, const T &b, const T &c)

{

return (b.x - a.x)*(c.y - a.y) - (b.y - a.y)*(c.x - a.x);

}

// 2D Signed Triangle Area

template<typename T>

inline float tri_area(const T &a, const T &b, const T &c)

{

return 0.5f*((b.x - a.x)*(c.y - a.y) - (c.x - a.x)*(b.y - a.y));

}

// 2D Signed Quad Area

template <typename T>

inline float quad_area(const T &a, const T &b, const T &c, const T &d)

{

return 0.5f*((c.x - a.x)*(d.y - b.y) - (d.x - b.x)*(c.y - a.y));

}

// Intersection test of two 2D line segments, returns intersect time along first segment, -1 if no intersection

template<typename T>

inline float segInt(const T &s1, const T &e1, const T &s2, const T &e2)

{

const float dx11 = e1.x - s1.x, dx22 = e2.x - s2.x;

const float dy11 = e1.y - s1.y, dy22 = e2.y - s2.y;

const float d = dy22*dx11 - dx22*dy11;

if (d == 0.f) return -1.f;

const float dx12 = s1.x - s2.x, dy12 = s1.y - s2.y;

const float a = dx22*dy12 - dy22*dx12, b = dx11*dy12 - dy11*dx12;

const float u1 = a / d;

if ((u1 < 0.f) || (u1 > 1.f)) return -1.f;

const float u2 = b / d;

return ((u2 < 0.f) || (u2 > 1.f)) ? -1.f : u1;

}

// Area of completely general quad

template<typename T>

inline float gQuad_area(const T &a, const T &b, const T &c, const T &d)

{

// test if any two opposite segments intersect,

// if so, return sum of area of the two triangles connected at intersection point,

// otherwise, do direct area computation which also works for concave quads.

float u = segInt(a, b, c, d);

if (u != -1.f) {

const T ip = a + u*(b - a);

return ABS(tri_area(a, ip, d)) + ABS(tri_area(b, ip, c));

}

u = segInt(a, d, b, c);

if (u != -1.f) {

const T ip = a + u*(d - a);

return ABS(tri_area(a, b, ip)) + ABS(tri_area(c, d, ip));

}

return ABS(quad_area(a, b, c, d));

}

// returns minimum of two elements

template<typename T>

inline T MIN2(T a, T b)

{

return (a < b) ? a : b;

}

// returns maximum of two elements

template<typename T>

inline T MAX2(T a, T b)

{

return (a > b) ? a : b;

}

// returns minimum of four elements

template<typename T>

inline T MIN4(T a, T b, T c, T d)

100

{

101

return MIN2(a, MIN2(b, MIN2(c, d)));

102

}

103

104

// returns maximum of four elements

105

template<typename T>

106

inline T MAX4(T a, T b, T c, T d)

107

{

108

return MAX2(a, MAX2(b, MAX2(c, d)));

109

}

110

111

// swaps contents of two objects

112

template<typename T>

113

inline void SWAP(T &a, T &b)

114

{

115

T t = a;

116

a = b;

117

b = t;

118

}

119

120

// clamp value to zero if negative

121

template<typename T>

122

inline T CLAMP0(T v)

123

{

124

return (v < (T)0) ? ((T)0) : v;

125

}

126

127

// clamp value in range [min, max]

128

template<typename T>

129

inline T CLAMP(T a, T min, T max)

130

{

131

return ((a < min) ? min : ((a > max) ? max : a));

132

}

133

134

// returns sign of x

135

template<typename T>

136

inline T SIGN(T x) { return (x==(T)0) ? 0 : ((x<(T)0) ? (T)-1 : (T)1); }

137

138

//----------------------------------------------------------------------------------

139

// General simple C based math routines mainly used by the shader virtual machine,

140

// most vector functions are also used for points/normals and colors.

141

// For the cases where the difference matters, explicit routines exist,

142

// eg. matrix vector & point multiply .etc.

143

// General format is func(result, var1, var2, ...)

144

// For all functions where operands are the same type as the result type,

145

// the result variable can be the same as either of the operands, eg. cross(v1, v1, v2)

146

//----------------------------------------------------------------------------------

147

// vectors

148

inline RtVoid addVVV(RtVector r, const RtVector v1, const RtVector v2)

149

{

150

r[0] = v1[0] + v2[0];

151

r[1] = v1[1] + v2[1];

152

r[2] = v1[2] + v2[2];

153

}

154

155

inline RtVoid subVVV(RtVector r, const RtVector v1, const RtVector v2)

156

{

157

r[0] = v1[0] - v2[0];

158

r[1] = v1[1] - v2[1];

159

r[2] = v1[2] - v2[2];

160

}

161

162

inline RtVoid mulVVV(RtVector r, const RtVector v1, const RtVector v2)

163

{

164

r[0] = v1[0] * v2[0];

165

r[1] = v1[1] * v2[1];

166

r[2] = v1[2] * v2[2];

167

}

168

169

inline RtVoid maddVVV(RtVector r, const RtVector v1, const RtVector v2)

170

{

171

r[0] += v1[0] * v2[0];

172

r[1] += v1[1] * v2[1];

173

r[2] += v1[2] * v2[2];

174

}

175

176

inline RtVoid msubVVV(RtVector r, const RtVector v1, const RtVector v2)

177

{

178

r[0] -= v1[0] * v2[0];

179

r[1] -= v1[1] * v2[1];

180

r[2] -= v1[2] * v2[2];

181

}

182

183

inline RtVoid divVVV(RtVector r, const RtVector v1, const RtVector v2)

184

{

185

if (v2[0] != 0.f) r[0] = v1[0] / v2[0];

186

if (v2[1] != 0.f) r[1] = v1[1] / v2[1];

187

if (v2[2] != 0.f) r[2] = v1[2] / v2[2];

188

}

189

190

inline RtVoid mulVVF(RtVector r, const RtVector v, const RtFloat f)

191

{

192

r[0] = v[0] * f;

193

r[1] = v[1] * f;

194

r[2] = v[2] * f;

195

}

196

197

inline RtVoid maddVVF(RtVector r, const RtVector v, const RtFloat f)

198

{

199

r[0] += v[0] * f;

200

r[1] += v[1] * f;

201

r[2] += v[2] * f;

202

}

203

204

inline RtVoid msubVVF(RtVector r, const RtVector v, const RtFloat f)

205

{

206

r[0] -= v[0] * f;

207

r[1] -= v[1] * f;

208

r[2] -= v[2] * f;

209

}

210

211

inline RtVoid divVVF(RtVector r, const RtVector v, RtFloat f)

212

{

213

if (f != 0.f) {

214

f = 1.f/f;

215

r[0] = v[0] * f;

216

r[1] = v[1] * f;

217

r[2] = v[2] * f;

218

}

219

}

220

221

inline RtVoid mulVV(RtVector r, const RtVector v)

222

{

223

r[0] *= v[0];

224

r[1] *= v[1];

225

r[2] *= v[2];

226

}

227

228

inline RtVoid vcross(RtVector r, const RtVector v1, const RtVector v2)

229

{

230

const RtFloat x = v1[1]*v2[2] - v1[2]*v2[1];

231

const RtFloat y = v1[2]*v2[0] - v1[0]*v2[2];

232

r[2] = v1[0]*v2[1] - v1[1]*v2[0];

233

r[0] = x;

234

r[1] = y;

235

}

236

237

inline RtVoid vdot(RtFloat &r, const RtVector v1, const RtVector v2)

238

{

239

r = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];

240

}

241

242

inline RtVoid vnormalize(RtVector r, const RtVector v)

243

{

244

RtFloat d = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];

245

if (d != 0.f) {

246

d = 1.f / sqrtf(d);

247

r[0] = v[0]*d;

248

r[1] = v[1]*d;

249

r[2] = v[2]*d;

250

}

251

}

252

253

// sets the length of v to L

254

inline RtVoid vsetlength(RtVector r, const RtVector v, float L)

255

{

256

RtFloat d = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];

257

if (d != 0.f) {

258

d = L / sqrtf(d);

259

r[0] = v[0]*d;

260

r[1] = v[1]*d;

261

r[2] = v[2]*d;

262

}

263

}

264

265

inline RtVoid vlength(RtFloat &r, const RtVector v)

266

{

267

r = sqrtf(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);

268

}

269

270

// step/smoothstep

271

inline RtVoid step(RtFloat &r, const RtFloat min, const RtFloat val)

272

{

273

r = (val < min) ? 0.f : 1.f;

274

}

275

276

inline RtVoid smoothstep(RtFloat &r, const RtFloat min, const RtFloat max, RtFloat val)

277

{

278

if (val < min)

279

r = 0.f;

280

else if (val > max)

281

r = 1.f;

282

else {

283

val = (val - min) / (max - min);

284

r = val*val*(3.f - 2.f*val);

285

}

286

}

287

288

// faceforward

289

inline RtVoid faceforward(RtVector r, const RtVector v1, const RtVector v2, const RtVector v3)

290

{

291

// assumes v1 = N, I = v2, v3 = Ng|Nref

292

const RtFloat d = ((-v2[0]*v3[0] - v2[1]*v3[1] - v2[2]*v3[2]) < 0.f) ? -1.f : 1.f;

293

r[0] = d*v1[0];

294

r[1] = d*v1[1];

295

r[2] = d*v1[2];

296

}

297

298

// float linear interpolation

299

inline RtVoid mixf(RtFloat &r, const RtFloat f0, const RtFloat f1, const RtFloat v)

300

{

301

r = f0 + v*(f1 - f0);

302

}

303

304

// color/point/vector/normal linear interpolation

305

inline RtVoid mixv(RtVector r, const RtVector c0, const RtVector c1, const RtFloat v)

306

{

307

r[0] = c0[0] + v*(c1[0] - c0[0]);

308

r[1] = c0[1] + v*(c1[1] - c0[1]);

309

r[2] = c0[2] + v*(c1[2] - c0[2]);

310

}

311

312

// conversion of degrees to radians

313

inline RtVoid radians(RtFloat &r, const RtFloat d)

314

{

315

r = d * ((RtFloat)M_PI/180.f);

316

}

317

318

// conversion of radians to degrees

319

inline RtVoid degrees(RtFloat &r, const RtFloat d)

320

{

321

r = d * (180.f/(RtFloat)M_PI);

322

}

323

324

// bilinear interpolation (other float[3] types can be used too, RtPoint/RtVector/etc)

325

inline void bilerp(RtColor r, RtFloat u, RtFloat v,

326

const RtColor c00, const RtColor c10, const RtColor c01, const RtColor c11)

327

{

328

const RtFloat w00=(1.f-u)*(1.f-v), w10=u*(1.f-v), w01=(1.f-u)*v, w11=u*v;

329

r[0] = w00*c00[0] + w10*c10[0] + w01*c01[0] + w11*c11[0];

330

r[1] = w00*c00[1] + w10*c10[1] + w01*c01[1] + w11*c11[1];

331

r[2] = w00*c00[2] + w10*c10[2] + w01*c01[2] + w11*c11[2];

332

}

333

334

// bilinear interpolation of single floats

335

inline void bilerpF(RtFloat &r, RtFloat u, RtFloat v,

336

RtFloat f00, RtFloat f10, RtFloat f01, RtFloat f11)

337

{

338

r = (1.f-v)*((1.f-u)*f00 + u*f10) + v*((1.f-u)*f01 + u*f11);

339

}

340

341

// reflect vector

342

inline void reflect(RtVector r, const RtVector I, const RtVector N)

343

{

344

const RtFloat ndi2 = 2.f*(N[0]*I[0] + N[1]*I[1] + N[2]*I[2]);

345

r[0] = I[0] - ndi2*N[0];

346

r[1] = I[1] - ndi2*N[1];

347

r[2] = I[2] - ndi2*N[2];

348

}

349

350

// refract vector

351

inline void refract(RtVector r, const RtVector I, const RtVector N, const RtFloat eta)

352

{

353

const RtFloat ndi = N[0]*I[0] + N[1]*I[1] + N[2]*I[2];

354

const RtFloat k = 1.f - eta*eta*(1.f - ndi*ndi);

355

if (k < 0.f)

356

r[0] = r[1] = r[2] = 0.f;

357

else {

358

const RtFloat f = eta*ndi + sqrtf(k);

359

r[0] = eta*I[0] - f*N[0];

360

r[1] = eta*I[1] - f*N[1];

361

r[2] = eta*I[2] - f*N[2];

362

}

363

}

364

365

void fresnel(const RtVector I, const RtVector N, const float eta, float& Kr, float& Kt, RtVector R = NULL, RtVector T = NULL);

366

367

//------------------------------------------------------------------------------

368

// matrices

369

370

RtVoid addMMM(RtMatrix r, const RtMatrix m1, const RtMatrix m2);

371

RtVoid subMMM(RtMatrix r, const RtMatrix m1, const RtMatrix m2);

372

RtVoid mulMMM(RtMatrix r, const RtMatrix m1, const RtMatrix m2);

373

RtBoolean invertMatrix(RtMatrix minv, const RtMatrix m2);

374

375

inline RtVoid transposeM(RtMatrix r)

376

{

377

SWAP(r[0][1], r[1][0]);

378

SWAP(r[0][2], r[2][0]);

379

SWAP(r[0][3], r[3][0]);

380

SWAP(r[2][1], r[1][2]);

381

SWAP(r[3][2], r[2][3]);

382

SWAP(r[3][1], r[1][3]);

383

}

384

385

inline RtVoid divMMM(RtMatrix r, const RtMatrix m1, const RtMatrix m2)

386

{

387

RtMatrix im2;

388

assert(invertMatrix(im2, m2));

389

mulMMM(r, m1, im2);

390

}

391

392

inline RtVoid mulVMV(RtVector r, const RtMatrix m, const RtVector v)

393

{

394

const RtFloat x = v[0], y = v[1];

395

r[0] = m[0][0]*x + m[0][1]*y + m[0][2]*v[2];

396

r[1] = m[1][0]*x + m[1][1]*y + m[1][2]*v[2];

397

r[2] = m[2][0]*x + m[2][1]*y + m[2][2]*v[2];

398

}

399

400

inline RtVoid mulNMN(RtNormal r, const RtMatrix m, const RtNormal n)

401

{

402

// normal, multiply by inverse of transpose (could be optimized, don't need to recalculate everytime)

403

RtMatrix m2;

404

assert(invertMatrix(m2, m));

405

const RtFloat x = n[0], y = n[1];

406

r[0] = m2[0][0]*x + m2[1][0]*y + m2[2][0]*n[2];

407

r[1] = m2[0][1]*x + m2[1][1]*y + m2[2][1]*n[2];

408

r[2] = m2[0][2]*x + m2[1][2]*y + m2[2][2]*n[2];

409

}

410

411

inline RtVoid mulPMP(RtPoint r, const RtMatrix m, const RtPoint p)

412

{

413

const RtFloat x = p[0], y = p[1], z = p[2];

414

r[0] = m[0][0]*x + m[0][1]*y + m[0][2]*z + m[0][3];

415

r[1] = m[1][0]*x + m[1][1]*y + m[1][2]*z + m[1][3];

416

r[2] = m[2][0]*x + m[2][1]*y + m[2][2]*z + m[2][3];

417

const RtFloat w = m[3][0]*x + m[3][1]*y + m[3][2]*z + m[3][3];

418

if (w != 1.f) divVVF(r, r, w);

419

}

420

421

inline RtVoid mulPMP4(RtHpoint r, const RtMatrix m, const RtHpoint p)

422

{

423

const RtFloat x = p[0], y = p[1], z = p[2];

424

r[0] = m[0][0]*x + m[0][1]*y + m[0][2]*z + m[0][3]*p[3];

425

r[1] = m[1][0]*x + m[1][1]*y + m[1][2]*z + m[1][3]*p[3];

426

r[2] = m[2][0]*x + m[2][1]*y + m[2][2]*z + m[2][3]*p[3];

427

r[3] = m[3][0]*x + m[3][1]*y + m[3][2]*z + m[3][3]*p[3];

428

}

429

430

//------------------------------------------------------------------------------

431

// Not sure what the purpose of the below floor/ceil/fmod functions was anymore,

432

// kind of superfluous, could just as well use the actual math funcs...

433

434

// returns highest integer <= x as float

435

inline float FLOORF(float x)

436

{

437

const float r = float(int(x));

438

return ((x >= 0.f) || (r == x)) ? r : (r - 1.f);

439

}

440

441

// returns lowest integer >= x as float

442

inline float CEILF(float x)

443

{

444

const float r = float(int(x));

445

return ((x <= 0.f) || (r == x)) ? r : (r + 1.f);

446

}

447

448

// returns highest integer <= x as integer

449

inline int FLOORI(float x)

450

{

451

if (x >= 0.f) return int(x);

452

const int r = int(x);

453

return (float(r) == x) ? r : (r - 1);

454

}

455

456

// returns lowest integer >= x as integer

457

inline int CEILI(float x)

458

{

459

if (x <= 0.f) return int(x);

460

const int r = int(x);

461

return (float(r) == x) ? r : (r + 1);

462

}

463

464

// returns float x modulo y, return value is in range (0, y)

465

inline float FMOD(float x, float y)

466

{

467

float z = x / y;

468

return y*(z - FLOORF(z));

469

}

470

471

// returns integer log2 of x (== highest set bit)

472

inline unsigned int ilog2(unsigned int x)

473

{

474

unsigned int b = 0;

475

while (x >>= 1) b++;

476

return b;

477

}

478

479

// as above, but always rounds upwards

480

inline unsigned int ilog2_roundup(unsigned int x)

481

{

482

unsigned int xx = x, b = 0;

483

while (xx >>= 1) b++;

484

return b + (unsigned int)((unsigned int)(1 << b) < x);

485

}

486

487

// direct bicubic evaluation funcs

488

// direct Bezier eval. at value t

489

inline Point3 Bezier(float t, const Point3 &p0, const Point3 &p1, const Point3 &p2, const Point3 &p3)

490

{

491

const float tm = 1.f - t;

492

return (3.f*tm*p1 + 3.f*t*p2)*tm*t + tm*tm*tm*p0 + t*t*t*p3;

493

}

494

495

// derivative of Bezier

496

inline void BezierDeriv(float t, const RtPoint p0, const RtPoint p1, const RtPoint p2, const RtPoint p3, RtPoint bd)

497

{

498

// for now the more complicated form here

499

bd[0] = (-p0[0] + 3.f*p1[0] + (-3.f*p2[0]) + p3[0])*3.f*t*t + (3.f*p0[0] + (-6.f*p1[0]) + 3.f*p2[0])*2.f*t + (-3.f*p0[0] + 3.f*p1[0]);

500

bd[1] = (-p0[1] + 3.f*p1[1] + (-3.f*p2[1]) + p3[1])*3.f*t*t + (3.f*p0[1] + (-6.f*p1[1]) + 3.f*p2[1])*2.f*t + (-3.f*p0[1] + 3.f*p1[1]);

501

bd[2] = (-p0[2] + 3.f*p1[2] + (-3.f*p2[2]) + p3[2])*3.f*t*t + (3.f*p0[2] + (-6.f*p1[2]) + 3.f*p2[2])*2.f*t + (-3.f*p0[2] + 3.f*p1[2]);

502

}

503

504

// direct evaluation using Power basis

505

inline Point3 Power(float t, const Point3 &p0, const Point3 &p1, const Point3 &p2, const Point3 &p3)

506

{

507

return ((p0*t + p1)*t + p2)*t + p3;

508

}

509

510

// direct bicubic curve eval. at value t, arbitrary basis b

511

inline void SplineF(RtFloat& r, const float t, const RtBasis &b, const RtFloat &f0, const RtFloat &f1, const RtFloat &f2, const RtFloat &f3)

512

{

513

r = (((b[0][0]*f0 + b[0][1]*f1 + b[0][2]*f2 + b[0][3]*f3) *t

514

+ (b[1][0]*f0 + b[1][1]*f1 + b[1][2]*f2 + b[1][3]*f3))*t

515

+ (b[2][0]*f0 + b[2][1]*f1 + b[2][2]*f2 + b[2][3]*f3))*t

516

+ (b[3][0]*f0 + b[3][1]*f1 + b[3][2]*f2 + b[3][3]*f3);

517

}

518

519

inline void SplineV(RtVector r, const float t, const RtBasis &b, const RtVector &v0, const RtVector &v1, const RtVector &v2, const RtVector &v3)

520

{

521

SplineF(r[0], t, b, v0[0], v1[0], v2[0], v3[0]);

522

SplineF(r[1], t, b, v0[1], v1[1], v2[1], v3[1]);

523

SplineF(r[2], t, b, v0[2], v1[2], v2[2], v3[2]);

524

}

525

526

// Bezier eval. using de Casteljau method

527

inline void deCasteljau_P(float t, const RtPoint* cp, RtPoint P)

528

{

529

// osculating plane -> tangent -> point

530

float d0, d1, d2;

531

#define CAS(a, b, c, d, r) \

532

d0 = a + t*(b - a); \

533

d1 = b + t*(c - b); \

534

d2 = c + t*(d - c); \

535

d0 += t*(d1 - d0); \

536

d1 += t*(d2 - d1); \

537

r = d0 + t*(d1 - d0);

538

CAS(cp[0][0], cp[1][0], cp[2][0], cp[3][0], P[0])

539

CAS(cp[0][1], cp[1][1], cp[2][1], cp[3][1], P[1])

540

CAS(cp[0][2], cp[1][2], cp[2][2], cp[3][2], P[2])

541

#undef CAS

542

}

543

544

// as above, but also returns v tangent in dPdv

545

inline void deCasteljau_P_dPdv(float t, const RtPoint* cp, RtPoint P, RtVector dPdv)

546

{

547

float d0, d1, d2;

548

#define CAS(a, b, c, d, r, rt) \

549

d0 = a + t*(b - a); \

550

d1 = b + t*(c - b); \

551

d2 = c + t*(d - c); \

552

d0 += t*(d1 - d0); \

553

d1 += t*(d2 - d1); \

554

r = d0 + t*(rt = d1 - d0); \

555

rt *= 3.f;

556

CAS(cp[0][0], cp[1][0], cp[2][0], cp[3][0], P[0], dPdv[0])

557

CAS(cp[0][1], cp[1][1], cp[2][1], cp[3][1], P[1], dPdv[1])

558

CAS(cp[0][2], cp[1][2], cp[2][2], cp[3][2], P[2], dPdv[2])

559

#undef CAS

560

}

561

562

// as above, but only returns u tangent vector instead

563

inline void deCasteljau_dPdu(float t, const RtPoint* cp, RtVector dPdu)

564

{

565

float d0, d1, d2;

566

#define CAS(a, b, c, d, r) \

567

d0 = a + t*(b - a); \

568

d1 = b + t*(c - b); \

569

d2 = c + t*(d - c); \

570

d0 += t*(d1 - d0); \

571

d1 += t*(d2 - d1); \

572

r = 3.f*(d1 - d0);

573

CAS(cp[0][0], cp[1][0], cp[2][0], cp[3][0], dPdu[0])

574

CAS(cp[0][1], cp[1][1], cp[2][1], cp[3][1], dPdu[1])

575

CAS(cp[0][2], cp[1][2], cp[2][2], cp[3][2], dPdu[2])

576

#undef CAS

577

}

578

579

// version of deCasteljau_P() which can used for c.pts not in array

580

inline void deCasteljau_P(float t, const RtPoint cp0, const RtPoint cp1,

581

const RtPoint cp2, const RtPoint cp3, RtPoint P)

582

{

583

float d0, d1, d2;

584

#define CAS(a, b, c, d, r) \

585

d0 = a + t*(b - a); \

586

d1 = b + t*(c - b); \

587

d2 = c + t*(d - c); \

588

d0 += t*(d1 - d0); \

589

d1 += t*(d2 - d1); \

590

r = d0 + t*(d1 - d0);

591

CAS(cp0[0], cp1[0], cp2[0], cp3[0], P[0])

592

CAS(cp0[1], cp1[1], cp2[1], cp3[1], P[1])

593

CAS(cp0[2], cp1[2], cp2[2], cp3[2], P[2])

594

#undef CAS

595

}

596

597

__END_QDRENDER

598

599

#endif //_MATHUTIL_H

Older »