1*cd705ea8SMichael Walle // SPDX-License-Identifier: GPL-2.0-only
2*cd705ea8SMichael Walle /*
3*cd705ea8SMichael Walle * Generic polynomial calculation using integer coefficients.
4*cd705ea8SMichael Walle *
5*cd705ea8SMichael Walle * Copyright (C) 2020 BAIKAL ELECTRONICS, JSC
6*cd705ea8SMichael Walle *
7*cd705ea8SMichael Walle * Authors:
8*cd705ea8SMichael Walle * Maxim Kaurkin <maxim.kaurkin@baikalelectronics.ru>
9*cd705ea8SMichael Walle * Serge Semin <Sergey.Semin@baikalelectronics.ru>
10*cd705ea8SMichael Walle *
11*cd705ea8SMichael Walle */
12*cd705ea8SMichael Walle
13*cd705ea8SMichael Walle #include <linux/kernel.h>
14*cd705ea8SMichael Walle #include <linux/module.h>
15*cd705ea8SMichael Walle #include <linux/polynomial.h>
16*cd705ea8SMichael Walle
17*cd705ea8SMichael Walle /*
18*cd705ea8SMichael Walle * Originally this was part of drivers/hwmon/bt1-pvt.c.
19*cd705ea8SMichael Walle * There the following conversion is used and should serve as an example here:
20*cd705ea8SMichael Walle *
21*cd705ea8SMichael Walle * The original translation formulae of the temperature (in degrees of Celsius)
22*cd705ea8SMichael Walle * to PVT data and vice-versa are following:
23*cd705ea8SMichael Walle *
24*cd705ea8SMichael Walle * N = 1.8322e-8*(T^4) + 2.343e-5*(T^3) + 8.7018e-3*(T^2) + 3.9269*(T^1) +
25*cd705ea8SMichael Walle * 1.7204e2
26*cd705ea8SMichael Walle * T = -1.6743e-11*(N^4) + 8.1542e-8*(N^3) + -1.8201e-4*(N^2) +
27*cd705ea8SMichael Walle * 3.1020e-1*(N^1) - 4.838e1
28*cd705ea8SMichael Walle *
29*cd705ea8SMichael Walle * where T = [-48.380, 147.438]C and N = [0, 1023].
30*cd705ea8SMichael Walle *
31*cd705ea8SMichael Walle * They must be accordingly altered to be suitable for the integer arithmetics.
32*cd705ea8SMichael Walle * The technique is called 'factor redistribution', which just makes sure the
33*cd705ea8SMichael Walle * multiplications and divisions are made so to have a result of the operations
34*cd705ea8SMichael Walle * within the integer numbers limit. In addition we need to translate the
35*cd705ea8SMichael Walle * formulae to accept millidegrees of Celsius. Here what they look like after
36*cd705ea8SMichael Walle * the alterations:
37*cd705ea8SMichael Walle *
38*cd705ea8SMichael Walle * N = (18322e-20*(T^4) + 2343e-13*(T^3) + 87018e-9*(T^2) + 39269e-3*T +
39*cd705ea8SMichael Walle * 17204e2) / 1e4
40*cd705ea8SMichael Walle * T = -16743e-12*(D^4) + 81542e-9*(D^3) - 182010e-6*(D^2) + 310200e-3*D -
41*cd705ea8SMichael Walle * 48380
42*cd705ea8SMichael Walle * where T = [-48380, 147438] mC and N = [0, 1023].
43*cd705ea8SMichael Walle *
44*cd705ea8SMichael Walle * static const struct polynomial poly_temp_to_N = {
45*cd705ea8SMichael Walle * .total_divider = 10000,
46*cd705ea8SMichael Walle * .terms = {
47*cd705ea8SMichael Walle * {4, 18322, 10000, 10000},
48*cd705ea8SMichael Walle * {3, 2343, 10000, 10},
49*cd705ea8SMichael Walle * {2, 87018, 10000, 10},
50*cd705ea8SMichael Walle * {1, 39269, 1000, 1},
51*cd705ea8SMichael Walle * {0, 1720400, 1, 1}
52*cd705ea8SMichael Walle * }
53*cd705ea8SMichael Walle * };
54*cd705ea8SMichael Walle *
55*cd705ea8SMichael Walle * static const struct polynomial poly_N_to_temp = {
56*cd705ea8SMichael Walle * .total_divider = 1,
57*cd705ea8SMichael Walle * .terms = {
58*cd705ea8SMichael Walle * {4, -16743, 1000, 1},
59*cd705ea8SMichael Walle * {3, 81542, 1000, 1},
60*cd705ea8SMichael Walle * {2, -182010, 1000, 1},
61*cd705ea8SMichael Walle * {1, 310200, 1000, 1},
62*cd705ea8SMichael Walle * {0, -48380, 1, 1}
63*cd705ea8SMichael Walle * }
64*cd705ea8SMichael Walle * };
65*cd705ea8SMichael Walle */
66*cd705ea8SMichael Walle
67*cd705ea8SMichael Walle /**
68*cd705ea8SMichael Walle * polynomial_calc - calculate a polynomial using integer arithmetic
69*cd705ea8SMichael Walle *
70*cd705ea8SMichael Walle * @poly: pointer to the descriptor of the polynomial
71*cd705ea8SMichael Walle * @data: input value of the polynimal
72*cd705ea8SMichael Walle *
73*cd705ea8SMichael Walle * Calculate the result of a polynomial using only integer arithmetic. For
74*cd705ea8SMichael Walle * this to work without too much loss of precision the coefficients has to
75*cd705ea8SMichael Walle * be altered. This is called factor redistribution.
76*cd705ea8SMichael Walle *
77*cd705ea8SMichael Walle * Returns the result of the polynomial calculation.
78*cd705ea8SMichael Walle */
polynomial_calc(const struct polynomial * poly,long data)79*cd705ea8SMichael Walle long polynomial_calc(const struct polynomial *poly, long data)
80*cd705ea8SMichael Walle {
81*cd705ea8SMichael Walle const struct polynomial_term *term = poly->terms;
82*cd705ea8SMichael Walle long total_divider = poly->total_divider ?: 1;
83*cd705ea8SMichael Walle long tmp, ret = 0;
84*cd705ea8SMichael Walle int deg;
85*cd705ea8SMichael Walle
86*cd705ea8SMichael Walle /*
87*cd705ea8SMichael Walle * Here is the polynomial calculation function, which performs the
88*cd705ea8SMichael Walle * redistributed terms calculations. It's pretty straightforward.
89*cd705ea8SMichael Walle * We walk over each degree term up to the free one, and perform
90*cd705ea8SMichael Walle * the redistributed multiplication of the term coefficient, its
91*cd705ea8SMichael Walle * divider (as for the rationale fraction representation), data
92*cd705ea8SMichael Walle * power and the rational fraction divider leftover. Then all of
93*cd705ea8SMichael Walle * this is collected in a total sum variable, which value is
94*cd705ea8SMichael Walle * normalized by the total divider before being returned.
95*cd705ea8SMichael Walle */
96*cd705ea8SMichael Walle do {
97*cd705ea8SMichael Walle tmp = term->coef;
98*cd705ea8SMichael Walle for (deg = 0; deg < term->deg; ++deg)
99*cd705ea8SMichael Walle tmp = mult_frac(tmp, data, term->divider);
100*cd705ea8SMichael Walle ret += tmp / term->divider_leftover;
101*cd705ea8SMichael Walle } while ((term++)->deg);
102*cd705ea8SMichael Walle
103*cd705ea8SMichael Walle return ret / total_divider;
104*cd705ea8SMichael Walle }
105*cd705ea8SMichael Walle EXPORT_SYMBOL_GPL(polynomial_calc);
106*cd705ea8SMichael Walle
107*cd705ea8SMichael Walle MODULE_DESCRIPTION("Generic polynomial calculations");
108*cd705ea8SMichael Walle MODULE_LICENSE("GPL");
109