1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
|
/*
* Copyright 2012 Google Inc.
*
* Use of this source code is governed by a BSD-style license that can be
* found in the LICENSE file.
*/
#include "DataTypes.h"
#include "Extrema.h"
static int validUnitDivide(double numer, double denom, double* ratio)
{
if (numer < 0) {
numer = -numer;
denom = -denom;
}
if (denom == 0 || numer == 0 || numer >= denom)
return 0;
double r = numer / denom;
if (r == 0) { // catch underflow if numer <<<< denom
return 0;
}
*ratio = r;
return 1;
}
/** From Numerical Recipes in C.
Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C])
x1 = Q / A
x2 = C / Q
*/
static int findUnitQuadRoots(double A, double B, double C, double roots[2])
{
if (A == 0)
return validUnitDivide(-C, B, roots);
double* r = roots;
double R = B*B - 4*A*C;
if (R < 0) { // complex roots
return 0;
}
R = sqrt(R);
double Q = (B < 0) ? -(B-R)/2 : -(B+R)/2;
r += validUnitDivide(Q, A, r);
r += validUnitDivide(C, Q, r);
if (r - roots == 2 && AlmostEqualUlps(roots[0], roots[1])) { // nearly-equal?
r -= 1; // skip the double root
}
return (int)(r - roots);
}
/** Cubic'(t) = At^2 + Bt + C, where
A = 3(-a + 3(b - c) + d)
B = 6(a - 2b + c)
C = 3(b - a)
Solve for t, keeping only those that fit between 0 < t < 1
*/
int findExtrema(double a, double b, double c, double d, double tValues[2])
{
// we divide A,B,C by 3 to simplify
double A = d - a + 3*(b - c);
double B = 2*(a - b - b + c);
double C = b - a;
return findUnitQuadRoots(A, B, C, tValues);
}
/** Quad'(t) = At + B, where
A = 2(a - 2b + c)
B = 2(b - a)
Solve for t, only if it fits between 0 < t < 1
*/
int findExtrema(double a, double b, double c, double tValue[1])
{
/* At + B == 0
t = -B / A
*/
return validUnitDivide(a - b, a - b - b + c, tValue);
}
|
