lib/polynomial.c

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