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
validUnitDivide(double numer,double denom,double * ratio)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 */
findUnitQuadRoots(double A,double B,double C,double roots[2])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 */
findExtrema(double a,double b,double c,double d,double tValues[2])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 */
findExtrema(double a,double b,double c,double tValue[1])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 }
82