symbolic4
integral.c
Go to the documentation of this file.
1 
2 /*
3 
4  Copyright (c) 2019 Hannes Eberhard
5 
6  Permission is hereby granted, free of charge, to any person obtaining a copy
7  of this software and associated documentation files (the "Software"), to deal
8  in the Software without restriction, including without limitation the rights
9  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10  copies of the Software, and to permit persons to whom the Software is
11  furnished to do so, subject to the following conditions:
12 
13  The above copyright notice and this permission notice shall be included in all
14  copies or substantial portions of the Software.
15 
16  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22  SOFTWARE.
23 
24  */
25 
26 #include "symbolic4.h"
27 
28 bool expression_is_risch_integrable(expression* source, expression* variable) {
29 
30  uint8_t i;
31 
32  for (i = 0; i < source->child_count; i++) {
33  if (!expression_is_risch_integrable(source->children[i], variable)) return false;
34  }
35 
36  if (source->identifier == EXPI_EXPONENTATION && source->children[0]->value.symbolic[0] != 'e' && count_occurrences(source->children[1], variable, true) != 0) {
37  return false;
38  } else if (source->identifier == EXPI_EXPONENTATION && count_occurrences(source->children[0], variable, true) != 0 && source->children[1]->identifier != EXPI_LITERAL) {
39  return false;
40  } else if (source->identifier == EXPI_EXPONENTATION && count_occurrences(source->children[0], variable, true) != 0 && source->children[1]->value.numeric.denominator != 1) {
41  return false;
42  } else if (source->identifier == EXPI_LOG) {
43  return false;
44  }
45 
46  return true;
47 
48 }
49 
50 void risch_get_extensions(expression* extensions, expression* source, expression* variable) {
51 
52  uint8_t i;
53  expression* extension;
54 
55  if (count_occurrences(source, variable, true) == 0) {
56  return;
57  }
58 
59  for (i = 0; i < source->child_count; i++) {
60  if (source->children[i] == NULL) continue;
61  risch_get_extensions(extensions, source->children[i], variable);
62  }
63 
64  if (source->identifier == EXPI_EXPONENTATION && source->children[0]->value.symbolic[0] == 'e') {
65  extension = new_expression(EXPT_STRUCTURE, EXPI_EXTENSION, 2,
66  new_symbol(EXPI_SYMBOL, "EE____"),
67  copy_expression(source->children[1]));
68  } else if (source->identifier == EXPI_LN) {
69  extension = new_expression(EXPT_STRUCTURE, EXPI_EXTENSION, 2,
70  new_symbol(EXPI_SYMBOL, "LE____"),
71  copy_expression(source->children[1]));
72  } else {
73  return;
74  }
75 
76  itoa(extension->children[0]->value.symbolic + 2, extensions->child_count);
77  extension->children[0]->value.symbolic[strlen(extension->children[0]->value.symbolic)] = '_';
78  append_child(extensions, extension);
79 
80  replace_expression(source, copy_expression(extension->children[0]));
81 
82 }
83 
84 uint8_t risch_determine_parts(expression** polynominal_part, expression** rational_part, const expression* source, const expression* variable, const expression* extensions) {
85 
86  uint8_t i;
87  expression* quotient;
88  expression* remainder;
89  expression* temp;
90 
91  if (expression_is_reziprocal(source)) {
92 
93  *polynominal_part = NULL;
94 
95  *rational_part = new_expression(EXPT_STRUCTURE, EXPI_LIST, 2,
96  new_literal(1, 1, 1),
97  copy_expression(source));
98 
99 
100  (*rational_part)->children[1]->children[1]->sign = 1;
101 
102  simplify(*rational_part, true);
103 
104  ERROR_CHECK(any_expression_to_sparse_polynomial((*rational_part)->children[0], variable));
105  ERROR_CHECK(any_expression_to_sparse_polynomial((*rational_part)->children[1], variable));
106 
107  } else if (source->identifier == EXPI_MULTIPLICATION) {
108 
109  temp = new_expression(EXPT_STRUCTURE, EXPI_LIST, 2,
110  new_expression(EXPT_OPERATION, EXPI_MULTIPLICATION, 0),
111  new_expression(EXPT_OPERATION, EXPI_MULTIPLICATION, 0));
112 
113  for (i = 0; i < source->child_count; i++) {
114  if (expression_is_reziprocal(source->children[i])) {
115  append_child(temp->children[1], copy_expression(source->children[i]));
116  temp->children[1]->children[temp->children[1]->child_count - 1]->children[1]->sign = 1;
117  } else {
118  append_child(temp->children[0], copy_expression(source->children[i]));
119  }
120  }
121 
122  simplify(temp, true);
123 
124  ERROR_CHECK(any_expression_to_sparse_polynomial(temp->children[0], variable));
125  ERROR_CHECK(validate_sparse_polynomial(temp->children[0], false, false, false));
126 
127  ERROR_CHECK(any_expression_to_sparse_polynomial(temp->children[1], variable));
128  ERROR_CHECK(validate_sparse_polynomial(temp->children[1], false, false, false));
129 
130  poly_div(&quotient, &remainder, temp->children[0], temp->children[1], -1);
131 
132  *polynominal_part = quotient;
133 
134  if (expressions_are_identical(remainder, new_literal(1, 0, 1), false)) {
135  *rational_part = NULL;
136  } else {
137  *rational_part = new_expression(EXPT_STRUCTURE, EXPI_LIST, 2,
138  remainder,
139  temp->children[1]);
140  }
141 
142  } else {
143 
144  *polynominal_part = copy_expression(source);
145  *rational_part = NULL;
146  ERROR_CHECK(any_expression_to_sparse_polynomial(*polynominal_part, variable));
147  ERROR_CHECK(validate_sparse_polynomial(*polynominal_part, false, false, false));
148 
149  }
150 
151  return RETS_SUCCESS;
152 
153 }
154 
155 void risch_integrate_polynominal_part(expression* source) {
156 
157  uint8_t i;
158  expression* exponent;
159  expression* result;
160 
161  if (source == NULL) return;
162 
163  result = new_expression(EXPT_STRUCTURE, EXPI_POLYNOMIAL_SPARSE, 0);
164 
165  for (i = 0; i < source->child_count; i++) {
166 
167  exponent = new_expression(EXPT_OPERATION, EXPI_ADDITION, 2,
168  copy_expression(source->children[i]->children[0]),
169  new_literal(1, 1, 1));
170 
171  simplify(exponent, true);
172 
173  append_child(result, new_expression(EXPT_STRUCTURE, EXPI_LIST, 3,
174  copy_expression(exponent),
175  new_expression(EXPT_OPERATION, EXPI_DIVISION, 2,
176  copy_expression(source->children[i]->children[1]),
177  copy_expression(exponent)),
178  copy_expression(source->children[i]->children[2])));
179 
180  free_expression(exponent, false);
181 
182  }
183 
184  simplify(result, true);
185  replace_expression(source, result);
186 
187 }
188 
189 void rothstein_trager_method(expression* source) {
190 
191  uint8_t i, j;
192  expression* a = copy_expression(source->children[0]);
193  uint8_t a_degree;
194  expression* b = copy_expression(source->children[1]);
195  expression* b_derivative;
196  expression* b_derivative_negative;
197  expression* a_coefficients = new_expression(EXPT_STRUCTURE, EXPI_LIST, 0);
198  expression* b_coefficients = new_expression(EXPT_STRUCTURE, EXPI_LIST, 0);
199  expression* z = new_symbol(EXPI_SYMBOL, "EZ");
200  expression* extentions = new_expression(EXPT_STRUCTURE, EXPI_LIST, 0);
201  expression* resultant_poly;
202  expression* resultant;
203  expression* primitive_part;
204  expression* primitive_part_factors;
205  expression* factor;
206  uint8_t factor_degree;
207  expression* equation;
208  expression* temp;
209  expression* gcd;
210  expression* result = new_expression(EXPT_OPERATION, EXPI_ADDITION, 0);
211 
212  derivative(&b_derivative, b, NULL, true);
213  b_derivative_negative = new_expression(EXPT_OPERATION, EXPI_MULTIPLICATION, 2,
214  copy_expression(b_derivative),
215  new_literal(-1, 1, 1));
216  simplify(b_derivative_negative, true);
217 
218  simplify(b, true);
219  any_expression_to_sparse_polynomial(b, NULL);
220  any_expression_to_sparse_polynomial(b_derivative, NULL);
221  any_expression_to_sparse_polynomial(b_derivative_negative, NULL);
222 
223  append_child(extentions, new_expression(EXPT_OPERATION, EXPI_EQUATION, 2,
224  copy_expression(z),
225  copy_expression(z)));
226 
227  resultant_poly = new_expression(EXPT_STRUCTURE, EXPI_POLYNOMIAL_DENSE, 2,
228  copy_expression(z),
229  new_expression(EXPT_STRUCTURE, EXPI_LIST, 2,
230  copy_expression(a),
231  copy_expression(b_derivative_negative)));
232 
233  any_expression_to_expression(resultant_poly->children[1]->children[0]);
234  any_expression_to_expression(resultant_poly->children[1]->children[1]);
235 
236  a_degree = max(a->children[0]->children[0]->value.numeric.numerator, b_derivative->children[0]->children[0]->value.numeric.numerator);
237 
238  any_expression_to_dense_polynomial(a, NULL);
239  any_expression_to_dense_polynomial(b, NULL);
240  any_expression_to_dense_polynomial(b_derivative_negative, NULL);
241 
242  for (i = 0; i < a_degree + 1; i++) {
243  append_child(a_coefficients, new_expression(EXPT_STRUCTURE, EXPI_POLYNOMIAL_DENSE, 2,
244  copy_expression(z),
245  new_expression(EXPT_STRUCTURE, EXPI_LIST, 2,
246  copy_expression(a->children[1]->children[i]),
247  copy_expression(b_derivative_negative->children[1]->children[i]))));
248  }
249 
250  for (i = 0; i < b->children[1]->child_count; i++) {
251  append_child(b_coefficients, new_expression(EXPT_STRUCTURE, EXPI_POLYNOMIAL_DENSE, 2,
252  copy_expression(z),
253  new_expression(EXPT_STRUCTURE, EXPI_LIST, 1,
254  copy_expression(b->children[1]->children[i]))));
255  }
256 
257  calculate_resultant(&resultant, a_coefficients, b_coefficients);
258  simplify(resultant, true);
259  primitive_part = copy_expression(resultant);
260 
261  if (primitive_part != NULL) {
262  make_monic(primitive_part);
263  }
264 
265  factor_square_free(&primitive_part_factors, primitive_part);
266 
267  for (i = 0; i < primitive_part_factors->child_count; i++) {
268 
269  factor = primitive_part_factors->children[0];
270 
271  if (literal_to_double(factor->children[0]) == 1) {
272  continue;
273  }
274 
275  factor_degree = factor->children[1]->value.numeric.numerator;
276 
277  if (factor_degree <= 2) {
278 
279  any_expression_to_expression_recursive(factor->children[0]);
280  any_expression_to_expression_recursive(a);
281  any_expression_to_expression_recursive(b);
282  any_expression_to_expression_recursive(b_derivative);
283 
284  equation = new_expression(EXPT_OPERATION, EXPI_EQUATION, 2,
285  factor->children[0],
286  new_literal(1, 0, 1));
287  solve(equation, NULL);
288  embed_in_list_if_necessary(equation);
289 
290  for (j = 0; j < equation->child_count; j++) {
291 
292  temp = new_expression(EXPT_OPERATION, EXPI_ADDITION, 2,
293  new_expression(EXPT_OPERATION, EXPI_MULTIPLICATION, 3,
294  new_literal(-1, 1, 1),
295  copy_expression(equation->children[j]->children[1]),
296  copy_expression(b_derivative)),
297  copy_expression(a));
298 
299  simplify(temp, true);
300 
301  any_expression_to_sparse_polynomial(temp, NULL);
302  poly_gcd(&gcd, temp, copy_expression(source->children[1]));
303  any_expression_to_expression_recursive(gcd);
304 
305  append_child(result, new_expression(EXPT_OPERATION, EXPI_MULTIPLICATION, 2,
306  new_expression(EXPT_FUNCTION, EXPI_LN, 1, copy_expression(gcd)),
307  copy_expression(equation->children[j]->children[1])));
308 
309 
310  }
311 
312  }
313 
314  }
315 
316  simplify(result, true);
317 
318  replace_expression(source, result);
319 
320 }
321 
322 uint8_t risch_integrate_rational_part(expression* source) {
323 
324  if (source == NULL) return RETS_UNCHANGED;
325 
326  if (poly_is_square_free(source->children[1])) {
327  rothstein_trager_method(source);
328  return RETS_CHANGED;
329  } else {
330  return RETS_UNCHANGED;
331  }
332 
333 }
334 
335 uint8_t risch_integrate(expression* source, expression* variable) {
336 
337  expression* polynominal_part;
338  expression* rational_part;
339  expression* extensions = new_expression(EXPT_STRUCTURE, EXPI_LIST, 0);
340  expression* last_extension;
341  expression* temp_source = copy_expression(source);
342 
343  if (!expression_is_risch_integrable(source, variable)) return RETS_UNCHANGED;
344 
345  risch_get_extensions(extensions, temp_source, variable);
346  last_extension = (extensions->child_count) ? copy_expression(extensions->children[extensions->child_count - 1]) : NULL;
347 
348  ERROR_CHECK(risch_determine_parts(&polynominal_part, &rational_part, temp_source, (last_extension) ? last_extension->children[0] : variable, extensions));
349 
350  if (extensions->child_count == 0) {
351  risch_integrate_polynominal_part(polynominal_part);
352  risch_integrate_rational_part(rational_part);
353  } else {
354  return RETS_UNCHANGED;
355  }
356 
357  if (polynominal_part == NULL) {
358  replace_expression(source, rational_part);
359  } else if (rational_part == NULL) {
360  replace_expression(source, polynominal_part);
361  } else {
362  replace_expression(source, new_expression(EXPT_OPERATION, EXPI_ADDITION, 2,
363  polynominal_part,
364  rational_part));
365  }
366 
367  simplify(source, true);
368 
369  return RETS_SUCCESS;
370 
371 }
372 
373 uint8_t antiderivative(expression** result, expression* source, expression* variable, bool persistent) {
374 
375  uint8_t i;
376  expression* temp_source = copy_expression(source);
377 
378  if (variable == NULL) {
379  variable = guess_symbol(source, "", 0);
380  if (variable == NULL) return set_error(ERRD_SYNTAX, ERRI_ARGUMENTS, "");
381  }
382 
383  if (count_occurrences(source, variable, true) == 0) {
384  *result = new_expression(EXPT_OPERATION, EXPI_MULTIPLICATION, 2,
385  copy_expression(source),
386  copy_expression(variable));
387  simplify(*result, true);
388  return RETS_SUCCESS;
389  }
390 
391  if (temp_source->identifier == EXPI_ADDITION) {
392  for (i = 0; i < temp_source->child_count; i++) {
393  ERROR_CHECK(risch_integrate(temp_source->children[i], variable));
394  }
395  } else {
396  ERROR_CHECK(risch_integrate(temp_source, variable));
397  }
398 
399  simplify(temp_source, true);
400  *result = temp_source;
401 
402  return RETS_SUCCESS;
403 
404 }
405 
406 uint8_t definite_integral(expression** result, expression* source, expression* variable, expression* lower_bound, expression* upper_bound) {
407 
408  expression* temp_source;
409  expression* upper_bound_value;
410  expression* lower_bound_value;
411 
412  if (variable == NULL) {
413  variable = guess_symbol(source, "", 0);
414  }
415 
416  antiderivative(&temp_source, source, variable, true);
417 
418  upper_bound_value = copy_expression(temp_source);
419  lower_bound_value = copy_expression(temp_source);
420 
421  replace_occurences(upper_bound_value, variable, upper_bound);
422  replace_occurences(lower_bound_value, variable, lower_bound);
423 
424  *result = new_expression(EXPT_OPERATION, EXPI_SUBTRACTION, 2,
425  upper_bound_value,
426  lower_bound_value);
427 
428  ERROR_CHECK(simplify(*result, true));
429 
430  return RETS_SUCCESS;
431 
432 }
433 
ERRI_ARGUMENTS
Definition: foundation.h:58
make_monic
void make_monic(expression *source)
Definition: polynomial.c:580
EXPI_LIST
Definition: expression.h:90
EXPT_FUNCTION
Definition: expression.h:35
numeric_value::denominator
uintmax_t denominator
Definition: expression.h:109
symbolic4.h
risch_integrate_polynominal_part
void risch_integrate_polynominal_part(expression *source)
Definition: integral.c:155
expression_is_risch_integrable
bool expression_is_risch_integrable(expression *source, expression *variable)
Definition: integral.c:28
EXPI_LN
Definition: expression.h:56
embed_in_list_if_necessary
void embed_in_list_if_necessary(expression *source)
Definition: expression.c:417
EXPI_EXPONENTATION
Definition: expression.h:53
expression::symbolic
char symbolic[10]
Definition: expression.h:128
EXPI_POLYNOMIAL_DENSE
Definition: expression.h:89
EXPI_SYMBOL
Definition: expression.h:45
derivative
uint8_t derivative(expression **result, expression *source, expression *variable, bool persistent)
Definition: derivative.c:209
replace_expression
void replace_expression(expression *a, expression *b)
Definition: expression.c:309
expression
Definition: expression.h:112
expression::identifier
expression_identifier identifier
Definition: expression.h:115
antiderivative
uint8_t antiderivative(expression **result, expression *source, expression *variable, bool persistent)
Definition: integral.c:373
literal_to_double
double literal_to_double(expression *source)
Definition: expression.c:908
replace_occurences
void replace_occurences(expression *source, const expression *child, const expression *replacement)
Definition: expression.c:688
any_expression_to_dense_polynomial
return_status any_expression_to_dense_polynomial(expression *source, const expression *variable)
Definition: polynomial.c:67
EXPT_STRUCTURE
Definition: expression.h:36
risch_get_extensions
void risch_get_extensions(expression *extensions, expression *source, expression *variable)
Definition: integral.c:50
ERRD_SYNTAX
Definition: foundation.h:36
validate_sparse_polynomial
return_status validate_sparse_polynomial(expression *source, bool allow_decimal_exponents, bool allow_negative_exponents, bool allow_arbitrary_base)
Definition: polynomial.c:82
expression::sign
int8_t sign
Definition: expression.h:117
ERROR_CHECK
#define ERROR_CHECK(F)
Check if the return status of a function is RETS_ERROR. If so, return RETS_ERROR.
Definition: foundation.h:31
EXPI_LOG
Definition: expression.h:57
new_symbol
expression * new_symbol(expression_identifier identifier, const char *value)
Allocates and initializes a new symbol/variable expression.
Definition: expression.c:252
expression::value
union expression::@0 value
solve
uint8_t solve(expression *source, expression *variable)
Definition: solve.c:337
RETS_UNCHANGED
Definition: foundation.h:69
RETS_SUCCESS
Definition: foundation.h:67
EXPI_EXTENSION
Definition: expression.h:92
risch_determine_parts
uint8_t risch_determine_parts(expression **polynominal_part, expression **rational_part, const expression *source, const expression *variable, const expression *extensions)
Definition: integral.c:84
poly_is_square_free
bool poly_is_square_free(expression *source)
Definition: polynomial.c:560
expressions_are_identical
bool expressions_are_identical(const expression *a, expression *b, bool persistent)
Definition: expression.c:467
RETS_CHANGED
Definition: foundation.h:68
itoa
void itoa(char *buffer, uintmax_t source)
Definition: foundation.c:216
free_expression
void free_expression(expression *source, bool persistent)
Definition: expression.c:315
EXPI_MULTIPLICATION
Definition: expression.h:51
set_error
uint8_t set_error(error_domain domain, error_identifier identifier, const char *body)
Sets current_error to the arguments provided.
Definition: foundation.c:152
EXPI_POLYNOMIAL_SPARSE
Definition: expression.h:88
count_occurrences
uint8_t count_occurrences(const expression *haystack, expression *needle, bool persistent)
Definition: expression.c:664
any_expression_to_sparse_polynomial
return_status any_expression_to_sparse_polynomial(expression *source, const expression *variable)
Definition: polynomial.c:53
EXPT_OPERATION
Definition: expression.h:34
calculate_resultant
void calculate_resultant(expression **result, expression *a, expression *b)
Definition: matrix.c:185
any_expression_to_expression
void any_expression_to_expression(expression *source)
Definition: polynomial.c:33
expression::numeric
numeric_value numeric
Definition: expression.h:129
max
uintmax_t max(uintmax_t a, uintmax_t b)
Definition: math_foundation.c:32
copy_expression
expression * copy_expression(const expression *source)
Returns a deep copy of an expression.
Definition: expression.c:286
risch_integrate
uint8_t risch_integrate(expression *source, expression *variable)
Definition: integral.c:335
EXPI_SUBTRACTION
Definition: expression.h:50
poly_gcd
uint8_t poly_gcd(expression **gcd, const expression *a, const expression *b)
Definition: polynomial.c:603
poly_div
uint8_t poly_div(expression **quotient, expression **remainder, const expression *a, const expression *b, int8_t degree)
Definition: polynomial.c:467
simplify
uint8_t simplify(expression *source, bool recursive)
Definition: simplify.c:1117
any_expression_to_expression_recursive
void any_expression_to_expression_recursive(expression *source)
Definition: polynomial.c:44
definite_integral
uint8_t definite_integral(expression **result, expression *source, expression *variable, expression *lower_bound, expression *upper_bound)
Definition: integral.c:406
append_child
void append_child(expression *parent, expression *child)
Appends a child to an expression.
Definition: expression.c:383
new_expression
expression * new_expression(expression_type type, expression_identifier identifier, uint8_t child_count,...)
Allocates and initializes a new expression with the arguments provided.
Definition: expression.c:183
EXPI_EQUATION
Definition: expression.h:48
rothstein_trager_method
void rothstein_trager_method(expression *source)
Definition: integral.c:189
guess_symbol
expression * guess_symbol(const expression *source, const char *custom_priorities, uint8_t rank)
Definition: expression.c:830
new_literal
expression * new_literal(int8_t sign, uintmax_t numerator, uintmax_t denominator)
Allocates and initializes a new literal expression.
Definition: expression.c:225
EXPI_ADDITION
Definition: expression.h:49
factor_square_free
uint8_t factor_square_free(expression **factors, const expression *source)
Definition: polynomial.c:680
risch_integrate_rational_part
uint8_t risch_integrate_rational_part(expression *source)
Definition: integral.c:322
expression::child_count
uint8_t child_count
Definition: expression.h:119
expression::children
struct expression ** children
Definition: expression.h:124
expression_is_reziprocal
bool expression_is_reziprocal(const expression *source)
Definition: expression.c:636
numeric_value::numerator
uintmax_t numerator
Definition: expression.h:108
EXPI_DIVISION
Definition: expression.h:52
EXPI_LITERAL
Definition: expression.h:44