Line data Source code
1 : /*
2 : * libpoporon - polynomial.c
3 : *
4 : * This file is part of libpoporon.
5 : *
6 : * Author: Brian Armstrong (original, libcorrect)
7 : * SPDX-License-Identifier: BSD-3-Clause
8 : */
9 :
10 : #include "internal/polynomial.h"
11 :
12 108 : extern polynomial_t polynomial_create(uint32_t order)
13 : {
14 : polynomial_t polynomial;
15 :
16 108 : polynomial.coeff = (uint16_t *)pcalloc(order + 1, sizeof(uint16_t));
17 108 : polynomial.order = order;
18 :
19 108 : return polynomial;
20 : }
21 :
22 108 : extern void polynomial_destroy(polynomial_t polynomial)
23 : {
24 108 : pfree(polynomial.coeff);
25 108 : }
26 :
27 429 : extern void polynomial_mul(poporon_gf_t *gf, polynomial_t l, polynomial_t r, polynomial_t res)
28 : {
29 : uint32_t i, j, j_limit;
30 :
31 429 : pmemset(res.coeff, 0, sizeof(uint16_t) * (res.order + 1));
32 1450 : for (i = 0; i <= l.order; i++) {
33 1021 : if (i > res.order) {
34 0 : continue;
35 : }
36 :
37 1021 : j_limit = (r.order > res.order - i) ? res.order - i : r.order;
38 :
39 17990 : for (j = 0; j <= j_limit; j++) {
40 16969 : res.coeff[i + j] = field_add(res.coeff[i + j], field_mul(gf, l.coeff[i], r.coeff[j]));
41 : }
42 : }
43 429 : }
44 :
45 7 : extern void polynomial_mod(poporon_gf_t *gf, polynomial_t dividend, polynomial_t divisor, polynomial_t mod)
46 : {
47 : uint32_t i, j, q_order;
48 : uint16_t divisor_leading, q_coeff;
49 :
50 7 : if (mod.order < dividend.order) {
51 0 : return;
52 : }
53 :
54 7 : pmemcpy(mod.coeff, dividend.coeff, sizeof(uint16_t) * (dividend.order + 1));
55 :
56 7 : divisor_leading = gf->log[divisor.coeff[divisor.order]];
57 :
58 1568 : for (i = dividend.order; i > 0; i--) {
59 1568 : if (i < divisor.order) {
60 7 : break;
61 : }
62 :
63 1561 : if (mod.coeff[i] == 0) {
64 1113 : continue;
65 : }
66 :
67 448 : q_order = i - divisor.order;
68 448 : q_coeff = field_div_log(gf, gf->log[mod.coeff[i]], divisor_leading);
69 :
70 15232 : for (j = 0; j <= divisor.order; j++) {
71 14784 : if (divisor.coeff[j] == 0) {
72 0 : continue;
73 : }
74 :
75 14784 : mod.coeff[j + q_order] =
76 14784 : field_add(mod.coeff[j + q_order], field_mul_log_element(gf, gf->log[divisor.coeff[j]], q_coeff));
77 : }
78 : }
79 : }
80 :
81 21 : extern void polynomial_formal_derivative(poporon_gf_t *gf, polynomial_t poly, polynomial_t der)
82 : {
83 : uint32_t i;
84 :
85 21 : pmemset(der.coeff, 0, sizeof(uint16_t) * (der.order + 1));
86 :
87 205 : for (i = 0; i <= der.order; i++) {
88 184 : der.coeff[i] = field_sum(poly.coeff[i + 1], i + 1);
89 : }
90 21 : }
91 :
92 0 : extern uint16_t polynomial_eval(poporon_gf_t *gf, polynomial_t poly, uint16_t val)
93 : {
94 : uint32_t i;
95 : uint16_t res, val_exponentiated, val_log;
96 :
97 0 : if (val == 0) {
98 0 : return poly.coeff[0];
99 : }
100 :
101 0 : res = 0;
102 0 : val_exponentiated = gf->log[1];
103 0 : val_log = gf->log[val];
104 :
105 0 : for (i = 0; i <= poly.order; i++) {
106 0 : if (poly.coeff[i] != 0) {
107 0 : res = field_add(res, field_mul_log_element(gf, gf->log[poly.coeff[i]], val_exponentiated));
108 : }
109 :
110 0 : val_exponentiated = field_mul_log(gf, val_exponentiated, val_log);
111 : }
112 :
113 0 : return res;
114 : }
115 :
116 1168 : extern uint16_t polynomial_eval_lut(poporon_gf_t *gf, polynomial_t poly, const uint16_t *val_exp)
117 : {
118 : uint32_t i;
119 : uint16_t res;
120 :
121 1168 : if (val_exp[0] == 0) {
122 0 : return poly.coeff[0];
123 : }
124 :
125 1168 : res = 0;
126 213310 : for (i = 0; i <= poly.order; i++) {
127 212142 : if (poly.coeff[i] != 0) {
128 79120 : res = field_add(res, field_mul_log_element(gf, gf->log[poly.coeff[i]], val_exp[i]));
129 : }
130 : }
131 :
132 1168 : return res;
133 : }
134 :
135 5632 : extern uint16_t polynomial_eval_log_lut(poporon_gf_t *gf, polynomial_t poly_log, const uint16_t *val_exp)
136 : {
137 : uint16_t res;
138 : uint32_t i;
139 :
140 5632 : if (val_exp[0] == 0) {
141 22 : if (poly_log.coeff[0] == 0) {
142 0 : return 0;
143 : }
144 :
145 22 : return gf->exp[poly_log.coeff[0]];
146 : }
147 :
148 5610 : res = 0;
149 53040 : for (i = 0; i <= poly_log.order; i++) {
150 47430 : if (poly_log.coeff[i] != 0) {
151 47175 : res = field_add(res, field_mul_log_element(gf, poly_log.coeff[i], val_exp[i]));
152 : }
153 : }
154 :
155 5610 : return res;
156 : }
157 :
158 1728 : extern void polynomial_build_exp_lut(poporon_gf_t *gf, uint16_t val, uint32_t order, uint16_t *val_exp)
159 : {
160 : uint32_t i;
161 : uint16_t val_exponentiated, val_log;
162 :
163 1728 : val_exponentiated = gf->log[1];
164 1728 : val_log = gf->log[val];
165 :
166 99840 : for (i = 0; i <= order; i++) {
167 98112 : if (val == 0) {
168 192 : val_exp[i] = 0;
169 : } else {
170 97920 : val_exp[i] = val_exponentiated;
171 97920 : val_exponentiated = field_mul_log(gf, val_exponentiated, val_log);
172 : }
173 : }
174 1728 : }
175 :
176 2 : extern polynomial_t polynomial_init_from_roots(poporon_gf_t *gf, uint32_t nroots, uint16_t *roots, polynomial_t poly,
177 : polynomial_t *scratch)
178 : {
179 : polynomial_t l, r[2];
180 : uint32_t i, rcoeffres, nextrcoeff;
181 : uint16_t l_coeff[2];
182 :
183 2 : l.order = 1;
184 2 : l.coeff = l_coeff;
185 :
186 2 : r[0] = scratch[0];
187 2 : r[1] = scratch[1];
188 2 : rcoeffres = 0;
189 :
190 2 : r[rcoeffres].coeff[1] = 1;
191 2 : r[rcoeffres].coeff[0] = roots[0];
192 2 : r[rcoeffres].order = 1;
193 :
194 2 : l.coeff[1] = 1;
195 :
196 36 : for (i = 1; i < nroots; i++) {
197 34 : l.coeff[0] = roots[i];
198 34 : nextrcoeff = rcoeffres;
199 34 : rcoeffres = (rcoeffres + 1) % 2;
200 34 : r[rcoeffres].order = i + 1;
201 34 : polynomial_mul(gf, l, r[nextrcoeff], r[rcoeffres]);
202 : }
203 :
204 2 : pmemcpy(poly.coeff, r[rcoeffres].coeff, (nroots + 1) * sizeof(uint16_t));
205 2 : poly.order = nroots;
206 :
207 2 : return poly;
208 : }
209 :
210 18 : extern polynomial_t polynomial_create_from_roots(poporon_gf_t *gf, uint32_t nroots, uint16_t *roots)
211 : {
212 : polynomial_t poly, l, r[2];
213 : uint32_t rcoeffres, nextrcoeff, i;
214 :
215 18 : poly = polynomial_create(nroots);
216 :
217 18 : l.order = 1;
218 18 : l.coeff = (uint16_t *)pcalloc(2, sizeof(uint16_t));
219 :
220 18 : r[0].coeff = (uint16_t *)pcalloc(nroots + 1, sizeof(uint16_t));
221 18 : r[0].order = 0;
222 18 : r[1].coeff = (uint16_t *)pcalloc(nroots + 1, sizeof(uint16_t));
223 18 : r[1].order = 0;
224 18 : rcoeffres = 0;
225 :
226 18 : r[rcoeffres].coeff[0] = roots[0];
227 18 : r[rcoeffres].coeff[1] = 1;
228 18 : r[rcoeffres].order = 1;
229 :
230 18 : l.coeff[1] = 1;
231 :
232 392 : for (i = 1; i < nroots; i++) {
233 374 : l.coeff[0] = roots[i];
234 374 : nextrcoeff = rcoeffres;
235 374 : rcoeffres = (rcoeffres + 1) % 2;
236 374 : r[rcoeffres].order = i + 1;
237 374 : polynomial_mul(gf, l, r[nextrcoeff], r[rcoeffres]);
238 : }
239 :
240 18 : pmemcpy(poly.coeff, r[rcoeffres].coeff, (nroots + 1) * sizeof(uint16_t));
241 18 : poly.order = nroots;
242 :
243 18 : pfree(l.coeff);
244 18 : pfree(r[0].coeff);
245 18 : pfree(r[1].coeff);
246 :
247 18 : return poly;
248 : }
|
