symbolic4
matrix.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 expression* new_matrix(uint8_t rows, uint8_t colums) {
29 
30  uint8_t i, j;
31  expression* result = new_expression(EXPT_STRUCTURE, EXPI_MATRIX, 0);
32 
33  for (i = 0; i < rows; i++) {
34  append_child(result, new_expression(EXPT_STRUCTURE, EXPI_LIST, 0));
35  for (j = 0; j < colums; j++) {
36  append_child(result->children[i], new_expression(EXPT_STRUCTURE, EXPI_LIST, 0));
37  }
38  }
39 
40  result->child_count = rows;
41 
42  return result;
43 
44 }
45 
46 expression* new_sub_matrix(expression* matrix, uint8_t rows, uint8_t colums) {
47 
48  uint8_t i, j;
49  uint8_t a;
50  uint8_t b;
51  expression* result = new_matrix(matrix->child_count - 1, matrix->children[0]->child_count - 1);
52 
53  for (i = 0; i < matrix->child_count - 1; i++) {
54  for (j = 0; j < matrix->children[0]->child_count -1; j++) {
55  a = i;
56  if (i >= rows) {
57  a = i + 1;
58  }
59  b = i;
60  if (j >= colums) {
61  b = j + 1;
62  }
63  result->children[i]->children[j] = copy_expression(matrix->children[a]->children[b]);
64  }
65  }
66 
67  return result;
68 
69 }
70 
71 void sylvester_matrix(expression** matrix, expression* a, expression* b) {
72 
73  uint8_t i, j;
74  expression* result;
75  uint8_t a_degree = a->child_count - 1;
76  uint8_t b_degree = b->child_count - 1;
77  uint8_t size = a_degree + b_degree;
78 
79  result = new_matrix(size, size);
80 
81  for (i = 0; i < a_degree + 1; i++) {
82  for (j = 0; j < b_degree; j++) {
83  result->children[j]->children[i + j] = copy_expression(a->children[a->child_count - 1 - i]);
84  }
85  }
86 
87  for (i = 0; i < b_degree + 1; i++) {
88  for (j = 0; j < a_degree; j++) {
89  result->children[j + b_degree]->children[i + j] = copy_expression(b->children[b->child_count - 1 - i]);
90  }
91  }
92 
93  *matrix = result;
94 
95 }
96 
97 void gauss_determinant(expression** determinant, const expression* matrix) {
98 
99  uint8_t i, j;
100  uint8_t factor = 1;
101  uint8_t max = 0;
102  expression* sub_matrix;
103  expression* f;
104  expression* new_row = new_expression(EXPT_STRUCTURE, EXPI_LIST, 0);
105  expression* new_determinant = smart_alloc(1, sizeof(expression));
106  expression* result;
107 
108  if (matrix->child_count == 1 && matrix->children[0]->child_count == 1) {
109  *determinant = copy_expression(matrix->children[0]->children[0]);
110  return;
111  }
112 
113  for (i = 0; i < matrix->child_count; i++) {
114  if (!expressions_are_equivalent(matrix->children[i]->children[0], new_literal(1, 0, 1), false)) {
115  max = i;
116  break;
117  }
118  }
119 
120  if (expressions_are_equivalent(matrix->children[max]->children[0], new_literal(1, 0, 1), false)) {
121 
122  *determinant = new_literal(1, 0, 1);
123  return;
124 
125  } else {
126 
127  swap_expressions(matrix->children[0], matrix->children[max]);
128 
129  if (max != 0) {
130  factor *= -1;
131  }
132 
133  sub_matrix = new_sub_matrix(copy_expression(matrix), 0, 0);
134 
135  for (i = 1; i < matrix->child_count; i++) {
136 
137  f = new_expression(EXPT_OPERATION, EXPI_DIVISION, 2,
138  new_expression(EXPT_OPERATION, EXPI_MULTIPLICATION, 2,
139  new_literal(-1, 1, 1),
140  copy_expression(matrix->children[i]->children[0])),
141  copy_expression(matrix->children[0]->children[0]));
142 
143  replace_null_with_zero(f);
144  any_expression_to_expression_recursive(f);
145 
146  simplify(f, true);
147 
148  if (expression_contains_division(f) != -1) {
149  make_monic(f->children[1]);
150  }
151 
152  new_row = new_expression(EXPT_STRUCTURE, EXPI_LIST, 0);
153 
154  for (j = 1; j < matrix->child_count; j++) {
155  append_child(new_row, new_expression(EXPT_OPERATION, EXPI_ADDITION, 2,
156  new_expression(EXPT_OPERATION, EXPI_MULTIPLICATION, 2,
157  copy_expression(f),
158  copy_expression(matrix->children[0]->children[j])),
159  copy_expression(sub_matrix->children[i - 1]->children[j - 1])));
160  }
161 
162  replace_null_with_zero(new_row);
163  any_expression_to_expression_recursive(new_row);
164  simplify(new_row, true);
165 
166  sub_matrix->children[i - 1] = copy_expression(new_row);
167 
168  }
169 
170  gauss_determinant(&new_determinant, sub_matrix);
171 
172  result = new_expression(EXPT_OPERATION, EXPI_MULTIPLICATION, 3,
173  copy_expression(matrix->children[0]->children[0]),
174  copy_expression(new_determinant),
175  new_literal(factor, 1, 1));
176 
177  simplify(result, true);
178 
179  }
180 
181  *determinant = result;
182 
183 }
184 
185 void calculate_resultant(expression** result, expression* a, expression* b) {
186 
187  expression* matrix = new_expression(EXPT_STRUCTURE, EXPI_MATRIX, 0);
188  expression* determinant = new_expression(EXPT_NULL, EXPI_NULL, 0);
189 
190  if (a->child_count == 1) {
191  if (expression_is_numerical(a->children[0])) {
192  if (expressions_are_identical(a->children[0], new_literal(1, 0, 1), false)) {
193  *result = new_literal(1, 0, 1);
194  return;
195  } else {
196  if (b->child_count == 0) {
197  if (expression_is_numerical(b->children[0])) {
198  if (!expressions_are_identical(b->children[0], new_literal(1, 0, 1), false)) {
199  *result = new_literal(1, 1, 1);
200  return;
201  }
202  }
203  }
204  }
205  } else {
206  if (expressions_are_identical(a->children[0], new_literal(1, 0, 1), false)) {
207  *result = new_literal(1, 0, 1);
208  return;
209  } else {
210  if (b->child_count == 0) {
211  if (expression_is_numerical(b->children[0])) {
212  if (!expressions_are_identical(b->children[0], new_literal(1, 0, 1), false)) {
213  *result = new_literal(1, 1, 1);
214  return;
215  }
216  }
217  }
218  }
219  }
220  }
221 
222  if (b->child_count == 1) {
223  if (expression_is_numerical(b->children[0])) {
224  *result = new_literal(1, 0, 1);
225  return;
226  }
227  }
228 
229  sylvester_matrix(&matrix, a, b);
230  gauss_determinant(&determinant, matrix);
231  simplify(determinant, true);
232  any_expression_to_dense_polynomial(determinant, new_symbol(EXPI_SYMBOL, "EZ"));
233 
234  *result = determinant;
235 
236 }
237 
238 
make_monic
void make_monic(expression *source)
Definition: polynomial.c:580
EXPI_LIST
Definition: expression.h:90
EXPI_NULL
Definition: expression.h:42
expression_contains_division
int8_t expression_contains_division(const expression *source)
Definition: expression.c:622
expression_is_numerical
bool expression_is_numerical(const expression *source)
Definition: expression.c:646
symbolic4.h
EXPI_SYMBOL
Definition: expression.h:45
expression
Definition: expression.h:112
new_sub_matrix
expression * new_sub_matrix(expression *matrix, uint8_t rows, uint8_t colums)
Definition: matrix.c:46
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
new_symbol
expression * new_symbol(expression_identifier identifier, const char *value)
Allocates and initializes a new symbol/variable expression.
Definition: expression.c:252
gauss_determinant
void gauss_determinant(expression **determinant, const expression *matrix)
Definition: matrix.c:97
sylvester_matrix
void sylvester_matrix(expression **matrix, expression *a, expression *b)
Definition: matrix.c:71
expressions_are_identical
bool expressions_are_identical(const expression *a, expression *b, bool persistent)
Definition: expression.c:467
EXPI_MULTIPLICATION
Definition: expression.h:51
swap_expressions
uint8_t swap_expressions(expression *a, expression *b)
Definition: expression.c:718
EXPT_OPERATION
Definition: expression.h:34
calculate_resultant
void calculate_resultant(expression **result, expression *a, expression *b)
Definition: matrix.c:185
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
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
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
replace_null_with_zero
void replace_null_with_zero(expression *source)
Definition: expression.c:703
EXPI_MATRIX
Definition: expression.h:91
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
new_matrix
expression * new_matrix(uint8_t rows, uint8_t colums)
Definition: matrix.c:28
EXPT_NULL
Definition: expression.h:32
expressions_are_equivalent
bool expressions_are_equivalent(const expression *a, expression *b, bool persistent)
Definition: expression.c:526
expression::child_count
uint8_t child_count
Definition: expression.h:119
expression::children
struct expression ** children
Definition: expression.h:124
EXPI_DIVISION
Definition: expression.h:52
smart_alloc
void * smart_alloc(uint8_t length, size_t size)
Allocates and keeps track of memory.
Definition: foundation.c:59