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