1/*
2 * Copyright 2018 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
8#ifndef SkCubicSolver_DEFINED
9#define SkCubicSolver_DEFINED
10
11#include "include/core/SkTypes.h"
12#include "include/private/SkFloatingPoint.h"
13
14//#define CUBICMAP_TRACK_MAX_ERROR
15
16namespace SK_OPTS_NS {
17
18 static float eval_poly(float t, float b) {
19 return b;
20 }
21
22 template <typename... Rest>
23 static float eval_poly(float t, float m, float b, Rest... rest) {
24 return eval_poly(t, sk_fmaf(m,t,b), rest...);
25 }
26
27 inline float cubic_solver(float A, float B, float C, float D) {
28 #ifdef CUBICMAP_TRACK_MAX_ERROR
29 static int max_iters = 0;
30 #endif
31
32 #ifdef SK_DEBUG
33 auto valid = [](float t) {
34 return t >= 0 && t <= 1;
35 };
36 #endif
37
38 auto guess_nice_cubic_root = [](float a, float b, float c, float d) {
39 return -d;
40 };
41 float t = guess_nice_cubic_root(A, B, C, D);
42
43 int iters = 0;
44 const int MAX_ITERS = 8;
45 for (; iters < MAX_ITERS; ++iters) {
46 SkASSERT(valid(t));
47 float f = eval_poly(t, A,B,C,D); // f = At^3 + Bt^2 + Ct + D
48 if (sk_float_abs(f) <= 0.00005f) {
49 break;
50 }
51 float fp = eval_poly(t, 3*A, 2*B, C); // f' = 3At^2 + 2Bt + C
52 float fpp = eval_poly(t, 3*A+3*A, 2*B); // f'' = 6At + 2B
53
54 float numer = 2 * fp * f;
55 float denom = sk_fmaf(2*fp, fp, -(f*fpp));
56
57 t -= numer / denom;
58 }
59
60 #ifdef CUBICMAP_TRACK_MAX_ERROR
61 if (max_iters < iters) {
62 max_iters = iters;
63 SkDebugf("max_iters %d\n", max_iters);
64 }
65 #endif
66 SkASSERT(valid(t));
67 return t;
68 }
69
70} // namespace SK_OPTS_NS
71#endif
72