1 /*
2 Qalculate (library)
3
4 Copyright (C) 2003-2007, 2008, 2016-2019 Hanna Knutsson (hanna.knutsson@protonmail.com)
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10 */
11
12 #include "support.h"
13
14 #include "MathStructure.h"
15 #include "Calculator.h"
16 #include "BuiltinFunctions.h"
17 #include "Number.h"
18 #include "Function.h"
19 #include "Variable.h"
20 #include "Unit.h"
21 #include "Prefix.h"
22
23 #include <map>
24 #include <algorithm>
25
26 #include "MathStructure-support.h"
27
28 using std::string;
29 using std::cout;
30 using std::vector;
31 using std::endl;
32
flattenVector(MathStructure & mstruct) const33 MathStructure &MathStructure::flattenVector(MathStructure &mstruct) const {
34 if(!isVector()) {
35 mstruct = *this;
36 return mstruct;
37 }
38 MathStructure mstruct2;
39 mstruct.clearVector();
40 for(size_t i = 0; i < SIZE; i++) {
41 if(CHILD(i).isVector()) {
42 CHILD(i).flattenVector(mstruct2);
43 for(size_t i2 = 0; i2 < mstruct2.size(); i2++) {
44 mstruct.addChild(mstruct2[i2]);
45 }
46 } else {
47 mstruct.addChild(CHILD(i));
48 }
49 }
50 return mstruct;
51 }
rankVector(bool ascending)52 bool MathStructure::rankVector(bool ascending) {
53 vector<int> ranked;
54 vector<bool> ranked_equals_prev;
55 bool b;
56 for(size_t index = 0; index < SIZE; index++) {
57 b = false;
58 for(size_t i = 0; i < ranked.size(); i++) {
59 if(CALCULATOR->aborted()) return false;
60 ComparisonResult cmp = CHILD(index).compare(CHILD(ranked[i]));
61 if(COMPARISON_NOT_FULLY_KNOWN(cmp)) {
62 CALCULATOR->error(true, _("Unsolvable comparison at element %s when trying to rank vector."), i2s(index).c_str(), NULL);
63 return false;
64 }
65 if((ascending && cmp == COMPARISON_RESULT_GREATER) || cmp == COMPARISON_RESULT_EQUAL || (!ascending && cmp == COMPARISON_RESULT_LESS)) {
66 if(cmp == COMPARISON_RESULT_EQUAL) {
67 ranked.insert(ranked.begin() + i + 1, index);
68 ranked_equals_prev.insert(ranked_equals_prev.begin() + i + 1, true);
69 } else {
70 ranked.insert(ranked.begin() + i, index);
71 ranked_equals_prev.insert(ranked_equals_prev.begin() + i, false);
72 }
73 b = true;
74 break;
75 }
76 }
77 if(!b) {
78 ranked.push_back(index);
79 ranked_equals_prev.push_back(false);
80 }
81 }
82 int n_rep = 0;
83 for(long int i = (long int) ranked.size() - 1; i >= 0; i--) {
84 if(CALCULATOR->aborted()) return false;
85 if(ranked_equals_prev[i]) {
86 n_rep++;
87 } else {
88 if(n_rep) {
89 MathStructure v(i + 1 + n_rep, 1L, 0L);
90 v += i + 1;
91 v *= MathStructure(1, 2, 0);
92 for(; n_rep >= 0; n_rep--) {
93 CHILD(ranked[i + n_rep]) = v;
94 }
95 } else {
96 CHILD(ranked[i]).set(i + 1, 1L, 0L);
97 }
98 n_rep = 0;
99 }
100 }
101 return true;
102 }
sortVector(bool ascending)103 bool MathStructure::sortVector(bool ascending) {
104 vector<size_t> ranked_mstructs;
105 bool b;
106 for(size_t index = 0; index < SIZE; index++) {
107 b = false;
108 for(size_t i = 0; i < ranked_mstructs.size(); i++) {
109 if(CALCULATOR->aborted()) return false;
110 ComparisonResult cmp = CHILD(index).compare(*v_subs[ranked_mstructs[i]]);
111 if(COMPARISON_MIGHT_BE_LESS_OR_GREATER(cmp)) {
112 CALCULATOR->error(true, _("Unsolvable comparison at element %s when trying to sort vector."), i2s(index).c_str(), NULL);
113 return false;
114 }
115 if((ascending && COMPARISON_IS_EQUAL_OR_GREATER(cmp)) || (!ascending && COMPARISON_IS_EQUAL_OR_LESS(cmp))) {
116 ranked_mstructs.insert(ranked_mstructs.begin() + i, v_order[index]);
117 b = true;
118 break;
119 }
120 }
121 if(!b) {
122 ranked_mstructs.push_back(v_order[index]);
123 }
124 }
125 v_order = ranked_mstructs;
126 return true;
127 }
getRange(int start,int end,MathStructure & mstruct) const128 MathStructure &MathStructure::getRange(int start, int end, MathStructure &mstruct) const {
129 if(!isVector()) {
130 if(start > 1) {
131 mstruct.clear();
132 return mstruct;
133 } else {
134 mstruct = *this;
135 return mstruct;
136 }
137 }
138 if(start < 1) start = 1;
139 else if(start > (long int) SIZE) {
140 mstruct.clear();
141 return mstruct;
142 }
143 if(end < (int) 1 || end > (long int) SIZE) end = SIZE;
144 else if(end < start) end = start;
145 mstruct.clearVector();
146 for(; start <= end; start++) {
147 mstruct.addChild(CHILD(start - 1));
148 }
149 return mstruct;
150 }
151
resizeVector(size_t i,const MathStructure & mfill)152 void MathStructure::resizeVector(size_t i, const MathStructure &mfill) {
153 if(i > SIZE) {
154 while(i > SIZE) {
155 APPEND(mfill);
156 }
157 } else if(i < SIZE) {
158 REDUCE(i)
159 }
160 }
161
rows() const162 size_t MathStructure::rows() const {
163 if(m_type != STRUCT_VECTOR || SIZE == 0 || (SIZE == 1 && (!CHILD(0).isVector() || CHILD(0).size() == 0))) return 0;
164 return SIZE;
165 }
columns() const166 size_t MathStructure::columns() const {
167 if(m_type != STRUCT_VECTOR || SIZE == 0 || !CHILD(0).isVector()) return 0;
168 return CHILD(0).size();
169 }
getElement(size_t row,size_t column) const170 const MathStructure *MathStructure::getElement(size_t row, size_t column) const {
171 if(row == 0 || column == 0 || row > rows() || column > columns()) return NULL;
172 if(CHILD(row - 1).size() < column) return NULL;
173 return &CHILD(row - 1)[column - 1];
174 }
getElement(size_t row,size_t column)175 MathStructure *MathStructure::getElement(size_t row, size_t column) {
176 if(row == 0 || column == 0 || row > rows() || column > columns()) return NULL;
177 if(CHILD(row - 1).size() < column) return NULL;
178 return &CHILD(row - 1)[column - 1];
179 }
getArea(size_t r1,size_t c1,size_t r2,size_t c2,MathStructure & mstruct) const180 MathStructure &MathStructure::getArea(size_t r1, size_t c1, size_t r2, size_t c2, MathStructure &mstruct) const {
181 size_t r = rows();
182 size_t c = columns();
183 if(r1 < 1) r1 = 1;
184 else if(r1 > r) r1 = r;
185 if(c1 < 1) c1 = 1;
186 else if(c1 > c) c1 = c;
187 if(r2 < 1 || r2 > r) r2 = r;
188 else if(r2 < r1) r2 = r1;
189 if(c2 < 1 || c2 > c) c2 = c;
190 else if(c2 < c1) c2 = c1;
191 mstruct.clearMatrix(); mstruct.resizeMatrix(r2 - r1 + 1, c2 - c1 + 1, m_undefined);
192 for(size_t index_r = r1; index_r <= r2; index_r++) {
193 for(size_t index_c = c1; index_c <= c2; index_c++) {
194 mstruct[index_r - r1][index_c - c1] = CHILD(index_r - 1)[index_c - 1];
195 }
196 }
197 return mstruct;
198 }
rowToVector(size_t r,MathStructure & mstruct) const199 MathStructure &MathStructure::rowToVector(size_t r, MathStructure &mstruct) const {
200 if(r > rows()) {
201 mstruct = m_undefined;
202 return mstruct;
203 }
204 if(r < 1) r = 1;
205 mstruct = CHILD(r - 1);
206 return mstruct;
207 }
columnToVector(size_t c,MathStructure & mstruct) const208 MathStructure &MathStructure::columnToVector(size_t c, MathStructure &mstruct) const {
209 if(c > columns()) {
210 mstruct = m_undefined;
211 return mstruct;
212 }
213 if(c < 1) c = 1;
214 mstruct.clearVector();
215 for(size_t i = 0; i < SIZE; i++) {
216 mstruct.addChild(CHILD(i)[c - 1]);
217 }
218 return mstruct;
219 }
matrixToVector(MathStructure & mstruct) const220 MathStructure &MathStructure::matrixToVector(MathStructure &mstruct) const {
221 if(!isVector()) {
222 mstruct = *this;
223 return mstruct;
224 }
225 mstruct.clearVector();
226 for(size_t i = 0; i < SIZE; i++) {
227 if(CHILD(i).isVector()) {
228 for(size_t i2 = 0; i2 < CHILD(i).size(); i2++) {
229 mstruct.addChild(CHILD(i)[i2]);
230 }
231 } else {
232 mstruct.addChild(CHILD(i));
233 }
234 }
235 return mstruct;
236 }
setElement(const MathStructure & mstruct,size_t row,size_t column)237 void MathStructure::setElement(const MathStructure &mstruct, size_t row, size_t column) {
238 if(row > rows() || column > columns() || row < 1 || column < 1) return;
239 CHILD(row - 1)[column - 1] = mstruct;
240 CHILD(row - 1).childUpdated(column);
241 CHILD_UPDATED(row - 1);
242 }
addRows(size_t r,const MathStructure & mfill)243 void MathStructure::addRows(size_t r, const MathStructure &mfill) {
244 if(r == 0) return;
245 size_t cols = columns();
246 MathStructure mstruct; mstruct.clearVector();
247 mstruct.resizeVector(cols, mfill);
248 for(size_t i = 0; i < r; i++) {
249 APPEND(mstruct);
250 }
251 }
addColumns(size_t c,const MathStructure & mfill)252 void MathStructure::addColumns(size_t c, const MathStructure &mfill) {
253 if(c == 0) return;
254 for(size_t i = 0; i < SIZE; i++) {
255 if(CHILD(i).isVector()) {
256 for(size_t i2 = 0; i2 < c; i2++) {
257 CHILD(i).addChild(mfill);
258 }
259 }
260 }
261 CHILDREN_UPDATED;
262 }
addRow(const MathStructure & mfill)263 void MathStructure::addRow(const MathStructure &mfill) {
264 addRows(1, mfill);
265 }
addColumn(const MathStructure & mfill)266 void MathStructure::addColumn(const MathStructure &mfill) {
267 addColumns(1, mfill);
268 }
resizeMatrix(size_t r,size_t c,const MathStructure & mfill)269 void MathStructure::resizeMatrix(size_t r, size_t c, const MathStructure &mfill) {
270 if(r > SIZE) {
271 addRows(r - SIZE, mfill);
272 } else if(r != SIZE) {
273 REDUCE(r);
274 }
275 size_t cols = columns();
276 if(c > cols) {
277 addColumns(c - cols, mfill);
278 } else if(c != cols) {
279 for(size_t i = 0; i < SIZE; i++) {
280 CHILD(i).resizeVector(c, mfill);
281 }
282 }
283 }
matrixIsSquare() const284 bool MathStructure::matrixIsSquare() const {
285 return rows() == columns();
286 }
287
isNumericMatrix() const288 bool MathStructure::isNumericMatrix() const {
289 if(!isMatrix()) return false;
290 for(size_t index_r = 0; index_r < SIZE; index_r++) {
291 for(size_t index_c = 0; index_c < CHILD(index_r).size(); index_c++) {
292 if(!CHILD(index_r)[index_c].isNumber() || CHILD(index_r)[index_c].isInfinity()) return false;
293 }
294 }
295 return true;
296 }
297
298 //from GiNaC
pivot(size_t ro,size_t co,bool symbolic)299 int MathStructure::pivot(size_t ro, size_t co, bool symbolic) {
300
301 size_t k = ro;
302 if(symbolic) {
303 while((k < SIZE) && (CHILD(k)[co].isZero())) ++k;
304 } else {
305 size_t kmax = k + 1;
306 Number mmax(CHILD(kmax)[co].number());
307 mmax.abs();
308 while(kmax < SIZE) {
309 if(CHILD(kmax)[co].number().isNegative()) {
310 Number ntmp(CHILD(kmax)[co].number());
311 ntmp.negate();
312 if(ntmp.isGreaterThan(mmax)) {
313 mmax = ntmp;
314 k = kmax;
315 }
316 } else if(CHILD(kmax)[co].number().isGreaterThan(mmax)) {
317 mmax = CHILD(kmax)[co].number();
318 k = kmax;
319 }
320 ++kmax;
321 }
322 if(!mmax.isZero()) k = kmax;
323 }
324 if(k == SIZE) return -1;
325 if(k == ro) return 0;
326
327 SWAP_CHILDREN(ro, k)
328
329 return k;
330
331 }
332
333
334 //from GiNaC
determinant_minor(const MathStructure & mtrx,MathStructure & mdet,const EvaluationOptions & eo)335 void determinant_minor(const MathStructure &mtrx, MathStructure &mdet, const EvaluationOptions &eo) {
336 size_t n = mtrx.size();
337 if(n == 1) {
338 mdet = mtrx[0][0];
339 return;
340 }
341 if(n == 2) {
342 mdet = mtrx[0][0];
343 mdet.calculateMultiply(mtrx[1][1], eo);
344 mdet.add(mtrx[1][0], true);
345 mdet[mdet.size() - 1].calculateMultiply(mtrx[0][1], eo);
346 mdet[mdet.size() - 1].calculateNegate(eo);
347 mdet.calculateAddLast(eo);
348 return;
349 }
350 if(n == 3) {
351 mdet = mtrx[0][0];
352 mdet.calculateMultiply(mtrx[1][1], eo);
353 mdet.calculateMultiply(mtrx[2][2], eo);
354 mdet.add(mtrx[0][0], true);
355 mdet[mdet.size() - 1].calculateMultiply(mtrx[1][2], eo);
356 mdet[mdet.size() - 1].calculateMultiply(mtrx[2][1], eo);
357 mdet[mdet.size() - 1].calculateNegate(eo);
358 mdet.calculateAddLast(eo);
359 mdet.add(mtrx[0][1], true);
360 mdet[mdet.size() - 1].calculateMultiply(mtrx[1][0], eo);
361 mdet[mdet.size() - 1].calculateMultiply(mtrx[2][2], eo);
362 mdet[mdet.size() - 1].calculateNegate(eo);
363 mdet.calculateAddLast(eo);
364 mdet.add(mtrx[0][2], true);
365 mdet[mdet.size() - 1].calculateMultiply(mtrx[1][0], eo);
366 mdet[mdet.size() - 1].calculateMultiply(mtrx[2][1], eo);
367 mdet.calculateAddLast(eo);
368 mdet.add(mtrx[0][1], true);
369 mdet[mdet.size() - 1].calculateMultiply(mtrx[1][2], eo);
370 mdet[mdet.size() - 1].calculateMultiply(mtrx[2][0], eo);
371 mdet.calculateAddLast(eo);
372 mdet.add(mtrx[0][2], true);
373 mdet[mdet.size() - 1].calculateMultiply(mtrx[1][1], eo);
374 mdet[mdet.size() - 1].calculateMultiply(mtrx[2][0], eo);
375 mdet[mdet.size() - 1].calculateNegate(eo);
376 mdet.calculateAddLast(eo);
377 return;
378 }
379
380 std::vector<size_t> Pkey;
381 Pkey.reserve(n);
382 std::vector<size_t> Mkey;
383 Mkey.reserve(n - 1);
384 typedef std::map<std::vector<size_t>, class MathStructure> Rmap;
385 typedef std::map<std::vector<size_t>, class MathStructure>::value_type Rmap_value;
386 Rmap A;
387 Rmap B;
388 for(size_t r = 0; r < n; ++r) {
389 Pkey.erase(Pkey.begin(), Pkey.end());
390 Pkey.push_back(r);
391 A.insert(Rmap_value(Pkey, mtrx[r][n - 1]));
392 }
393 for(long int c = n - 2; c >= 0; --c) {
394 Pkey.erase(Pkey.begin(), Pkey.end());
395 Mkey.erase(Mkey.begin(), Mkey.end());
396 for(size_t i = 0; i < n - c; ++i) Pkey.push_back(i);
397 size_t fc = 0;
398 do {
399 mdet.clear();
400 for(size_t r = 0; r < n - c; ++r) {
401 if (mtrx[Pkey[r]][c].isZero()) continue;
402 Mkey.erase(Mkey.begin(), Mkey.end());
403 for(size_t i = 0; i < n - c; ++i) {
404 if(i != r) Mkey.push_back(Pkey[i]);
405 }
406 mdet.add(mtrx[Pkey[r]][c], true);
407 mdet[mdet.size() - 1].calculateMultiply(A[Mkey], eo);
408 if(r % 2) mdet[mdet.size() - 1].calculateNegate(eo);
409 mdet.calculateAddLast(eo);
410 }
411 if(!mdet.isZero()) B.insert(Rmap_value(Pkey, mdet));
412 for(fc = n - c; fc > 0; --fc) {
413 ++Pkey[fc-1];
414 if(Pkey[fc - 1] < fc + c) break;
415 }
416 if(fc < n - c && fc > 0) {
417 for(size_t j = fc; j < n - c; ++j) Pkey[j] = Pkey[j - 1] + 1;
418 }
419 } while(fc);
420 A = B;
421 B.clear();
422 }
423 return;
424 }
425
426 //from GiNaC
gaussianElimination(const EvaluationOptions & eo,bool det)427 int MathStructure::gaussianElimination(const EvaluationOptions &eo, bool det) {
428
429 if(!isMatrix()) return 0;
430 bool isnumeric = isNumericMatrix();
431
432 size_t m = rows();
433 size_t n = columns();
434 int sign = 1;
435
436 size_t r0 = 0;
437 for(size_t c0 = 0; c0 < n && r0 < m - 1; ++c0) {
438 int indx = pivot(r0, c0, true);
439 if(indx == -1) {
440 sign = 0;
441 if(det) return 0;
442 }
443 if(indx >= 0) {
444 if(indx > 0) sign = -sign;
445 for(size_t r2 = r0 + 1; r2 < m; ++r2) {
446 if(!CHILD(r2)[c0].isZero()) {
447 if(isnumeric) {
448 Number piv(CHILD(r2)[c0].number());
449 piv /= CHILD(r0)[c0].number();
450 for(size_t c = c0 + 1; c < n; ++c) {
451 CHILD(r2)[c].number() -= piv * CHILD(r0)[c].number();
452 }
453 } else {
454 MathStructure piv(CHILD(r2)[c0]);
455 piv.calculateDivide(CHILD(r0)[c0], eo);
456 for(size_t c = c0 + 1; c < n; ++c) {
457 CHILD(r2)[c].add(piv, true);
458 CHILD(r2)[c][CHILD(r2)[c].size() - 1].calculateMultiply(CHILD(r0)[c], eo);
459 CHILD(r2)[c][CHILD(r2)[c].size() - 1].calculateNegate(eo);
460 CHILD(r2)[c].calculateAddLast(eo);
461 }
462 }
463 }
464 for(size_t c = r0; c <= c0; ++c) CHILD(r2)[c].clear();
465 }
466 if(det) {
467 for(size_t c = r0 + 1; c < n; ++c) CHILD(r0)[c].clear();
468 }
469 ++r0;
470 }
471 }
472 for(size_t r = r0 + 1; r < m; ++r) {
473 for(size_t c = 0; c < n; ++c) CHILD(r)[c].clear();
474 }
475
476 return sign;
477 }
478
479 //from GiNaC
480 template <class It>
permutation_sign(It first,It last)481 int permutation_sign(It first, It last)
482 {
483 if (first == last)
484 return 0;
485 --last;
486 if (first == last)
487 return 0;
488 It flag = first;
489 int sign = 1;
490
491 do {
492 It i = last, other = last;
493 --other;
494 bool swapped = false;
495 while (i != first) {
496 if (*i < *other) {
497 std::iter_swap(other, i);
498 flag = other;
499 swapped = true;
500 sign = -sign;
501 } else if (!(*other < *i))
502 return 0;
503 --i; --other;
504 }
505 if (!swapped)
506 return sign;
507 ++flag;
508 if (flag == last)
509 return sign;
510 first = flag;
511 i = first; other = first;
512 ++other;
513 swapped = false;
514 while (i != last) {
515 if (*other < *i) {
516 std::iter_swap(i, other);
517 flag = other;
518 swapped = true;
519 sign = -sign;
520 } else if (!(*i < *other))
521 return 0;
522 ++i; ++other;
523 }
524 if (!swapped)
525 return sign;
526 last = flag;
527 --last;
528 } while (first != last);
529
530 return sign;
531 }
532
533 //from GiNaC
determinant(MathStructure & mstruct,const EvaluationOptions & eo) const534 MathStructure &MathStructure::determinant(MathStructure &mstruct, const EvaluationOptions &eo) const {
535
536 if(!matrixIsSquare()) {
537 CALCULATOR->error(true, _("The determinant can only be calculated for square matrices."), NULL);
538 mstruct = m_undefined;
539 return mstruct;
540 }
541
542 if(SIZE == 1) {
543 mstruct = CHILD(0)[0];
544 } else if(isNumericMatrix()) {
545
546 mstruct.set(1, 1, 0);
547 MathStructure mtmp(*this);
548 int sign = mtmp.gaussianElimination(eo, true);
549 for(size_t d = 0; d < SIZE; ++d) {
550 mstruct.number() *= mtmp[d][d].number();
551 }
552 mstruct.number() *= sign;
553
554 } else {
555
556 typedef std::pair<size_t, size_t> sizet_pair;
557 std::vector<sizet_pair> c_zeros;
558 for(size_t c = 0; c < CHILD(0).size(); ++c) {
559 size_t acc = 0;
560 for(size_t r = 0; r < SIZE; ++r) {
561 if(CHILD(r)[c].isZero()) ++acc;
562 }
563 c_zeros.push_back(sizet_pair(acc, c));
564 }
565
566 std::sort(c_zeros.begin(), c_zeros.end());
567 std::vector<size_t> pre_sort;
568 for(std::vector<sizet_pair>::const_iterator i = c_zeros.begin(); i != c_zeros.end(); ++i) {
569 pre_sort.push_back(i->second);
570 }
571 std::vector<size_t> pre_sort_test(pre_sort);
572
573 int sign = permutation_sign(pre_sort_test.begin(), pre_sort_test.end());
574
575 MathStructure result;
576 result.clearMatrix();
577 result.resizeMatrix(SIZE, CHILD(0).size(), m_zero);
578
579 size_t c = 0;
580 for(std::vector<size_t>::const_iterator i = pre_sort.begin(); i != pre_sort.end(); ++i,++c) {
581 for(size_t r = 0; r < SIZE; ++r) result[r][c] = CHILD(r)[(*i)];
582 }
583 mstruct.clear();
584
585 determinant_minor(result, mstruct, eo);
586
587 if(sign != 1) {
588 mstruct.calculateMultiply(sign, eo);
589 }
590
591 }
592
593 mstruct.mergePrecision(*this);
594 return mstruct;
595
596 }
597
permanent(MathStructure & mstruct,const EvaluationOptions & eo) const598 MathStructure &MathStructure::permanent(MathStructure &mstruct, const EvaluationOptions &eo) const {
599 if(!matrixIsSquare()) {
600 CALCULATOR->error(true, _("The permanent can only be calculated for square matrices."), NULL);
601 mstruct = m_undefined;
602 return mstruct;
603 }
604 if(b_approx) mstruct.setApproximate();
605 mstruct.setPrecision(i_precision);
606 if(SIZE == 1) {
607 if(CHILD(0).size() >= 1) {
608 mstruct = CHILD(0)[0];
609 }
610 } else if(SIZE == 2) {
611 mstruct = CHILD(0)[0];
612 if(IS_REAL(mstruct) && IS_REAL(CHILD(1)[1])) {
613 mstruct.number() *= CHILD(1)[1].number();
614 } else {
615 mstruct.calculateMultiply(CHILD(1)[1], eo);
616 }
617 if(IS_REAL(mstruct) && IS_REAL(CHILD(1)[0]) && IS_REAL(CHILD(0)[1])) {
618 mstruct.number() += CHILD(1)[0].number() * CHILD(0)[1].number();
619 } else {
620 MathStructure mtmp = CHILD(1)[0];
621 mtmp.calculateMultiply(CHILD(0)[1], eo);
622 mstruct.calculateAdd(mtmp, eo);
623 }
624 } else {
625 MathStructure mtrx;
626 mtrx.clearMatrix();
627 mtrx.resizeMatrix(SIZE - 1, CHILD(0).size() - 1, m_undefined);
628 for(size_t index_c = 0; index_c < CHILD(0).size(); index_c++) {
629 for(size_t index_r2 = 1; index_r2 < SIZE; index_r2++) {
630 for(size_t index_c2 = 0; index_c2 < CHILD(index_r2).size(); index_c2++) {
631 if(index_c2 > index_c) {
632 mtrx.setElement(CHILD(index_r2)[index_c2], index_r2, index_c2);
633 } else if(index_c2 < index_c) {
634 mtrx.setElement(CHILD(index_r2)[index_c2], index_r2, index_c2 + 1);
635 }
636 }
637 }
638 MathStructure mdet;
639 mtrx.permanent(mdet, eo);
640 if(IS_REAL(mdet) && IS_REAL(CHILD(0)[index_c])) {
641 mdet.number() *= CHILD(0)[index_c].number();
642 } else {
643 mdet.calculateMultiply(CHILD(0)[index_c], eo);
644 }
645 if(IS_REAL(mdet) && IS_REAL(mstruct)) {
646 mstruct.number() += mdet.number();
647 } else {
648 mstruct.calculateAdd(mdet, eo);
649 }
650 }
651 }
652 return mstruct;
653 }
setToIdentityMatrix(size_t n)654 void MathStructure::setToIdentityMatrix(size_t n) {
655 clearMatrix();
656 resizeMatrix(n, n, m_zero);
657 for(size_t i = 0; i < n; i++) {
658 CHILD(i)[i] = m_one;
659 }
660 }
getIdentityMatrix(MathStructure & mstruct) const661 MathStructure &MathStructure::getIdentityMatrix(MathStructure &mstruct) const {
662 mstruct.setToIdentityMatrix(columns());
663 return mstruct;
664 }
665
666 //modified algorithm from eigenmath
invertMatrix(const EvaluationOptions & eo)667 bool MathStructure::invertMatrix(const EvaluationOptions &eo) {
668
669 if(!matrixIsSquare()) return false;
670
671 if(isNumericMatrix()) {
672
673 int d, i, j, n = SIZE;
674
675 MathStructure idstruct;
676 Number mtmp;
677 idstruct.setToIdentityMatrix(n);
678 MathStructure mtrx(*this);
679
680 for(d = 0; d < n; d++) {
681 if(mtrx[d][d].isZero()) {
682 for (i = d + 1; i < n; i++) {
683 if(!mtrx[i][d].isZero()) break;
684 }
685
686 if(i == n) {
687 CALCULATOR->error(true, _("Inverse of singular matrix."), NULL);
688 return false;
689 }
690
691 mtrx[i].ref();
692 mtrx[d].ref();
693 MathStructure *mptr = &mtrx[d];
694 mtrx.setChild_nocopy(&mtrx[i], d + 1);
695 mtrx.setChild_nocopy(mptr, i + 1);
696
697 idstruct[i].ref();
698 idstruct[d].ref();
699 mptr = &idstruct[d];
700 idstruct.setChild_nocopy(&idstruct[i], d + 1);
701 idstruct.setChild_nocopy(mptr, i + 1);
702
703 }
704
705 mtmp = mtrx[d][d].number();
706 mtmp.recip();
707
708 for(j = 0; j < n; j++) {
709
710 if(j > d) {
711 mtrx[d][j].number() *= mtmp;
712 }
713
714 idstruct[d][j].number() *= mtmp;
715
716 }
717
718 for(i = 0; i < n; i++) {
719
720 if(i == d) continue;
721
722 mtmp = mtrx[i][d].number();
723 mtmp.negate();
724
725 for(j = 0; j < n; j++) {
726
727 if(j > d) {
728 mtrx[i][j].number() += mtrx[d][j].number() * mtmp;
729 }
730
731 idstruct[i][j].number() += idstruct[d][j].number() * mtmp;
732
733 }
734 }
735 }
736 set_nocopy(idstruct);
737 MERGE_APPROX_AND_PREC(idstruct)
738 } else {
739 MathStructure *mstruct = new MathStructure();
740 determinant(*mstruct, eo);
741 mstruct->calculateInverse(eo);
742 adjointMatrix(eo);
743 multiply_nocopy(mstruct, true);
744 calculateMultiplyLast(eo);
745 }
746 return true;
747 }
748
adjointMatrix(const EvaluationOptions & eo)749 bool MathStructure::adjointMatrix(const EvaluationOptions &eo) {
750 if(!matrixIsSquare()) return false;
751 if(SIZE == 1) {CHILD(0)[0].set(1, 1, 0); return true;}
752 MathStructure msave(*this);
753 for(size_t index_r = 0; index_r < SIZE; index_r++) {
754 for(size_t index_c = 0; index_c < CHILD(0).size(); index_c++) {
755 msave.cofactor(index_r + 1, index_c + 1, CHILD(index_r)[index_c], eo);
756 }
757 }
758 transposeMatrix();
759 return true;
760 }
transposeMatrix()761 bool MathStructure::transposeMatrix() {
762 MathStructure msave(*this);
763 resizeMatrix(CHILD(0).size(), SIZE, m_undefined);
764 for(size_t index_r = 0; index_r < SIZE; index_r++) {
765 for(size_t index_c = 0; index_c < CHILD(0).size(); index_c++) {
766 CHILD(index_r)[index_c] = msave[index_c][index_r];
767 }
768 }
769 return true;
770 }
cofactor(size_t r,size_t c,MathStructure & mstruct,const EvaluationOptions & eo) const771 MathStructure &MathStructure::cofactor(size_t r, size_t c, MathStructure &mstruct, const EvaluationOptions &eo) const {
772 if(r < 1) r = 1;
773 if(c < 1) c = 1;
774 if(SIZE == 1 || r > SIZE || c > CHILD(0).size()) {
775 mstruct = m_undefined;
776 return mstruct;
777 }
778 r--; c--;
779 mstruct.clearMatrix();
780 mstruct.resizeMatrix(SIZE - 1, CHILD(0).size() - 1, m_undefined);
781 for(size_t index_r = 0; index_r < SIZE; index_r++) {
782 if(index_r != r) {
783 for(size_t index_c = 0; index_c < CHILD(0).size(); index_c++) {
784 if(index_c > c) {
785 if(index_r > r) {
786 mstruct[index_r - 1][index_c - 1] = CHILD(index_r)[index_c];
787 } else {
788 mstruct[index_r][index_c - 1] = CHILD(index_r)[index_c];
789 }
790 } else if(index_c < c) {
791 if(index_r > r) {
792 mstruct[index_r - 1][index_c] = CHILD(index_r)[index_c];
793 } else {
794 mstruct[index_r][index_c] = CHILD(index_r)[index_c];
795 }
796 }
797 }
798 }
799 }
800 MathStructure mstruct2;
801 mstruct = mstruct.determinant(mstruct2, eo);
802 if((r + c) % 2 == 1) {
803 mstruct.calculateNegate(eo);
804 }
805 return mstruct;
806 }
807
calculate_userfunctions(MathStructure & m,const MathStructure & x_mstruct,const EvaluationOptions & eo)808 bool calculate_userfunctions(MathStructure &m, const MathStructure &x_mstruct, const EvaluationOptions &eo) {
809 bool b_ret = false;
810 for(size_t i = 0; i < m.size(); i++) {
811 if(calculate_userfunctions(m[i], x_mstruct, eo)) {
812 m.childUpdated(i + 1);
813 b_ret = true;
814 }
815 }
816 if(m.isFunction()) {
817 if(!m.contains(x_mstruct, true)) {
818 m.calculateFunctions(eo);
819 b_ret = true;
820 } else if(m.function()->subtype() == SUBTYPE_USER_FUNCTION && m.function()->condition().empty()) {
821 bool b = true;
822 for(size_t i = 0; i < ((UserFunction*) m.function())->countSubfunctions(); i++) {
823 if(((UserFunction*) m.function())->subfunctionPrecalculated(i + 1)) {
824 b = false;
825 break;
826 }
827 }
828 for(size_t i = 0; b && i < m.size(); i++) {
829 Argument *arg = m.function()->getArgumentDefinition(i + 1);
830 if(arg && arg->tests() && (arg->type() != ARGUMENT_TYPE_FREE || !arg->getCustomCondition().empty() || arg->rationalPolynomial() || arg->zeroForbidden() || (arg->handlesVector() && m[i].isVector())) && m[i].contains(x_mstruct, true)) {
831 b = false;
832 break;
833 }
834 }
835 if(b && m.calculateFunctions(eo, false)) {
836 calculate_userfunctions(m, x_mstruct, eo);
837 b_ret = true;
838 }
839 }
840 }
841 return b_ret;
842 }
generateVector(MathStructure x_mstruct,const MathStructure & min,const MathStructure & max,int steps,MathStructure * x_vector,const EvaluationOptions & eo) const843 MathStructure MathStructure::generateVector(MathStructure x_mstruct, const MathStructure &min, const MathStructure &max, int steps, MathStructure *x_vector, const EvaluationOptions &eo) const {
844 if(steps < 1) {
845 steps = 1;
846 }
847 MathStructure x_value(min);
848 MathStructure y_value;
849 MathStructure y_vector;
850 y_vector.clearVector();
851 if(steps > 1000000) {
852 CALCULATOR->error(true, _("Too many data points"), NULL);
853 return y_vector;
854 }
855 MathStructure step(max);
856 step.calculateSubtract(min, eo);
857 step.calculateDivide(steps - 1, eo);
858 if(!step.isNumber() || step.number().isNegative()) {
859 CALCULATOR->error(true, _("The selected min and max do not result in a positive, finite number of data points"), NULL);
860 return y_vector;
861 }
862 y_vector.resizeVector(steps, m_zero);
863 if(x_vector) x_vector->resizeVector(steps, m_zero);
864 MathStructure mthis(*this);
865 mthis.unformat();
866 calculate_userfunctions(mthis, x_mstruct, eo);
867 for(int i = 0; i < steps; i++) {
868 if(x_vector) {
869 (*x_vector)[i] = x_value;
870 }
871 y_value = mthis;
872 y_value.replace(x_mstruct, x_value);
873 y_value.eval(eo);
874 y_vector[i] = y_value;
875 if(x_value.isNumber()) x_value.number().add(step.number());
876 else x_value.calculateAdd(step, eo);
877 if(CALCULATOR->aborted()) {
878 y_vector.resizeVector(i, m_zero);
879 if(x_vector) x_vector->resizeVector(i, m_zero);
880 return y_vector;
881 }
882 }
883 return y_vector;
884 }
generateVector(MathStructure x_mstruct,const MathStructure & min,const MathStructure & max,const MathStructure & step,MathStructure * x_vector,const EvaluationOptions & eo) const885 MathStructure MathStructure::generateVector(MathStructure x_mstruct, const MathStructure &min, const MathStructure &max, const MathStructure &step, MathStructure *x_vector, const EvaluationOptions &eo) const {
886 MathStructure x_value(min);
887 MathStructure y_value;
888 MathStructure y_vector;
889 y_vector.clearVector();
890 if(min != max) {
891 MathStructure mtest(max);
892 mtest.calculateSubtract(min, eo);
893 if(!step.isZero()) mtest.calculateDivide(step, eo);
894 if(step.isZero() || !mtest.isNumber() || mtest.number().isNegative()) {
895 CALCULATOR->error(true, _("The selected min, max and step size do not result in a positive, finite number of data points"), NULL);
896 return y_vector;
897 } else if(mtest.number().isGreaterThan(1000000)) {
898 CALCULATOR->error(true, _("Too many data points"), NULL);
899 return y_vector;
900 }
901 mtest.number().round();
902 unsigned int steps = mtest.number().uintValue();
903 y_vector.resizeVector(steps, m_zero);
904 if(x_vector) x_vector->resizeVector(steps, m_zero);
905 }
906 MathStructure mthis(*this);
907 mthis.unformat();
908 calculate_userfunctions(mthis, x_mstruct, eo);
909 ComparisonResult cr = max.compare(x_value);
910 size_t i = 0;
911 while(COMPARISON_IS_EQUAL_OR_LESS(cr)) {
912 if(x_vector) {
913 if(i >= x_vector->size()) x_vector->addChild(x_value);
914 else (*x_vector)[i] = x_value;
915 }
916 y_value = mthis;
917 y_value.replace(x_mstruct, x_value);
918 y_value.eval(eo);
919 if(i >= y_vector.size()) y_vector.addChild(y_value);
920 else y_vector[i] = y_value;
921 if(x_value.isNumber()) x_value.number().add(step.number());
922 else x_value.calculateAdd(step, eo);
923 cr = max.compare(x_value);
924 if(CALCULATOR->aborted()) {
925 y_vector.resizeVector(i, m_zero);
926 if(x_vector) x_vector->resizeVector(i, m_zero);
927 return y_vector;
928 }
929 i++;
930 }
931 y_vector.resizeVector(i, m_zero);
932 if(x_vector) x_vector->resizeVector(i, m_zero);
933 return y_vector;
934 }
generateVector(MathStructure x_mstruct,const MathStructure & x_vector,const EvaluationOptions & eo) const935 MathStructure MathStructure::generateVector(MathStructure x_mstruct, const MathStructure &x_vector, const EvaluationOptions &eo) const {
936 MathStructure y_value;
937 MathStructure y_vector;
938 y_vector.clearVector();
939 MathStructure mthis(*this);
940 mthis.unformat();
941 calculate_userfunctions(mthis, x_mstruct, eo);
942 for(size_t i = 1; i <= x_vector.countChildren(); i++) {
943 if(CALCULATOR->aborted()) {
944 y_vector.clearVector();
945 return y_vector;
946 }
947 y_value = mthis;
948 y_value.replace(x_mstruct, x_vector.getChild(i));
949 y_value.eval(eo);
950 y_vector.addChild(y_value);
951 }
952 return y_vector;
953 }
954
955