ECF 1.5
Algebra.h
1/*
2Algebra.h
3Written by Walter O. Krawec
4
5Copyright (c) 2016 Walter O. Krawec
6
7Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
8
9The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
10
11THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12*/
13
14#ifndef _ALGEBRA_H_
15#define _ALGEBRA_H_
16
17#include <math.h>
18#include <stdlib.h>
19#include <iostream>
20#include <string>
21#include <vector>
22#include <fstream>
23#include <complex.h>
24
25#define Complex double _Complex
26
27namespace algebra
28{
29
31 {
32 public:
33 static Complex one() {Complex out; out = 1.0 + 0*I; return out; }
34 static Complex zero() {Complex out; out = 0 + 0*I; return out;}
35
36 static void sq(Complex& v)
37 {
38 Complex v2 = creal(v) - cimag(v)*I;
39 v = mult(v, v2);
40 }
41 static Complex construct(double theta)
42 {
43 Complex out = cos(theta) + sin(theta)*I;
44 return out;
45 }
46
47 static Complex power(Complex& base, int e)
48 {
49 return cpow(base, e);
50 Complex a;
51 double real = log( sqrt( creal(base)*creal(base)* + cimag(base)*cimag(base)));
52 double imag = atan2(cimag(base), creal(base));
53
54 real *= e;
55 imag *= e;
56
57 a = real + imag*I;
58
59 Complex result = construct(cimag(a));
60
61
62 Complex r2 = creal(result)*exp(real) + cimag(result)*real*I;
63
64 return r2;
65 }
66
67 static Complex add(Complex& a, Complex& b)
68 {
69 return a+b;
70 }
71 static Complex subtract(Complex& a, Complex& b)
72 {
73 return a-b;
74 }
75
76 static Complex mult(Complex& a, Complex& b)
77 {
78 return a*b;
79 }
80
81 static Complex invert(Complex& a)
82 {
83 Complex out;
84 double d = creal(a)*creal(a) + cimag(a)*cimag(a);
85
86 out = creal(a) / d - cimag(a)/d *I;
87
88 return out;
89 }
90 };
91
92 class matrix
93 {
94 public:
95 matrix()
96 {
97 elements = NULL;
98 noCleanup = false;
99 r = 0;
100 c = 0;
101 n = 0;
102 }
103 matrix(const char* _name)
104 {
105 elements = NULL;
106 noCleanup = false;
107 r = 0;
108 c = 0;
109 n = 0;
110 varName = std::string(_name);
111 }
112 matrix(const matrix& m)
113 {
114 elements = NULL;
115 r = m.r;
116 c = m.c;
117 n = r*c;
118
119 if(r == 0 || c == 0 || !m.elements)
120 return;
121
122 elements = new Complex[n];
123
124 noCleanup = false;
125
126 for(int i=0; i<n; ++i)
127 elements[i] = m.elements[i];
128 }
129 ~matrix() {cleanup();}
130
131 void saveToFile(const std::string& filename);
132 void saveToFile(std::ofstream& f);
133
134 bool loadFromFile(const std::string& filename);
135 bool loadFromFile(std::ifstream& f);
136
137 matrix& operator=(const matrix& m)
138 {
139 if(!elements || r != m.r || c != m.c)
140 {
141 cleanup();
142 r = m.r;
143 c = m.c;
144 n = r*c;
145
146 if(r == 0 || c == 0 || !m.elements)
147 {
148 throw std::string("operator= error.\n");
149 return *this;
150 }
151
152 elements = new Complex[n];
153 }
154
155 if(!m.elements)
156 {
157 throw std::string("operator= error 2.\n");
158 return *this;
159 }
160
161 for(int i=0; i<n; ++i)
162 elements[i] = m.elements[i];
163
164 noCleanup = false;
165
166 return *this;
167 }
168
169 void cleanup()
170 {
171 if(elements && !noCleanup)
172 {
173 delete[] elements;
174 }
175 elements = NULL;
176 r = 0;
177 c = 0;
178 n = 0;
179 }
180
181 void zero()
182 {
183 if(!elements)
184 return;
185
186 for(int i=0; i<n; ++i)
187 {
188 elements[i] = 0.0+0.0*I;
189 }
190 }
191
192 void ident()
193 {
194 if(!elements)
195 return;
196
197 zero();
198 int dim = getRows();
199 if(getCols() < getRows())
200 dim = getCols();
201
202 Complex one = 1.0+0.0*I;
203
204 for(int i=0; i<dim; ++i)
205 this->setDirect(i,i,one);
206 }
207
208 Complex trace()
209 {
210 Complex c = 0.0 + 0.0*I;
211 for(int i=0; i<getRows(); ++i)
212 {
213 c = c + (*this)(i,i);
214 }
215 return c;
216 }
217
218 void clear(Complex value)
219 {
220 if(!elements)
221 return;
222 for(int i=0; i<n; ++i)
223 elements[i] = value;
224 }
225
226 bool create(int _r, int _c, Complex* _elements)
227 {
228 cleanup();
229 r = _r;
230 c = _c;
231 n = r*c;
232
233 elements = _elements;
234 noCleanup=true;
235 return true;
236 }
237 bool create(int _r, int _c)
238 {
239 if(_r <= 0 || _c <= 0)
240 return false;
241 if(elements && r == _r && c == _c)
242 return true;
243
244 cleanup();
245 r = _r;
246 c = _c;
247 n = r*c;
248
249 elements = new Complex[r*c];
250
251 noCleanup = false;
252
253 return true;
254 }
255
256 bool create(int _r)
257 {
258 return create(_r, _r);
259 }
260
261 Complex& operator() (int row, int col)
262 {
263 int _r = row;
264 int _c = col;
265 if(!elements || _r >= r || _c >= c || _r < 0 || _c < 0)
266 {
267 std::cout << "Error 5 (" << elements << ", " << _r << ", " << r << ", " << _c << ", " << c << ")\n";
268 std::cout << "Name = " << varName << "\n";
269 throw(std::string("Algebra Error 5.1"));
270 //throw(5);
271 //return double();
272 }
273
274 return elements[_c*r + _r];
275 }
276
277 void random_unitary(const std::vector <double>& phi_v, const std::vector <double>& psi_v, const std::vector <double>& chi_v, double alpha);
278 bool eigenValues(std::vector <double>& values);
279 double safeLog(double x)
280 {
281 if(x < 0.0000001)
282 return 0;
283 return log(x)/log(2.0);
284 }
285 double entropy()
286 {
287 std::vector<double> ev;
288 double sum = 0;
289
290 eigenValues(ev);
291
292 for(int i=0; i < ev.size(); ++i)
293 {
294 sum = sum - ev[i] * safeLog(ev[i]);
295 }
296
297 return sum;
298 }
299
300
301 Complex getDirect(int row, int col)
302 {
303 int _r = row;
304 int _c = col;
305 if(!elements || _r >= r || _c >= c || _r < 0 || _c < 0)
306 {
307 std::cout << "Error 5 (" << elements << ", " << _r << ", " << r << ", " << _c << ", " << c << ")\n";
308 throw(std::string("Algebra Error 5.2"));
309 //throw(5);
310 //return double();
311 }
312
313 return elements[_c*r + _r];
314 }
315
316 void setDirect(int _n, Complex data)
317 {
318 if(!elements || _n >= n)
319 {
320// std::cout << "Error 5 (" << elements << ", " << _r << ", " << r << ", " << _c << ", " << c << ")\n";
321 throw(std::string("Algebra Error 5.3"));
322 //throw(5);
323 //return double();
324 }
325
326 elements[_n] = data;
327 }
328
329 void setDirect(int row, int col, Complex data)
330 {
331 int _r = row;
332 int _c = col;
333 if(!elements || _r >= r || _c >= c || _r < 0 || _c < 0)
334 {
335 std::cout << "Error 5 (" << elements << ", " << _r << ", " << r << ", " << _c << ", " << c << ")\n";
336 throw(std::string("Algebra Error 5.4"));
337 //throw(5);
338 //return double();
339 }
340
341 elements[_c*r + _r] = data;
342 }
343
344 void setDirect(int row, int col, double re, double im)
345 {
346 int _r = row;
347 int _c = col;
348 if(!elements || _r >= r || _c >= c || _r < 0 || _c < 0)
349 {
350 std::cout << "Error 5 (" << elements << ", " << _r << ", " << r << ", " << _c << ", " << c << ")\n";
351 throw(std::string("Algebra Error 5.4"));
352 //throw(5);
353 //return double();
354 }
355
356 elements[_c*r + _r] = re + im*I;
357 }
358 Complex get(int _n)
359 {
360 if(!elements || _n >= n)
361 {
362 std::cout << "Error 2 (" << elements << ", " << _n << ", " << n << ")\n";
363 return 0;
364 }
365
366 return elements[_n];
367 }
368
369 Complex* getElements()
370 {
371 return elements;
372 }
373
374 void print(std::ostream& f, bool printAll=false)
375 {
376 if(!elements)
377 return;
378 int rStart = 0, cStart=0;
379 // if(r > 10)
380 // rStart = r-10;
381 // if(c > 10)
382 // cStart = c-10;
383
384 int rows = alg_min(10, r);
385 int cols = alg_min(10, c);
386
387 if(printAll)
388 {
389 rows = r;
390 cols = c;
391 }
392 for(int _r=rStart; _r < rows; ++_r)
393 {
394 for(int _c=cStart; _c<cols; ++_c)
395 {
396 f << creal((*this)(_r,_c)) << "+" << cimag((*this)(_r,_c)) << "i\t";
397 }
398 f << "\n";
399 }
400
401 f << "\n";
402 }
403
404 void multiplyScalar(Complex& c)
405 {
406 if(!elements)
407 return;
408
409 for(int i=0; i<n; ++i)
410 {
411 elements[i] = ComplexHelper::mult(elements[i], c);
412 }
413 }
414
415 void sqTranspose()
416 {
417 if(!elements || c != r)
418 {
419 throw(std::string("Algebra error 1.2"));
420 return;
421 }
422
423 Complex temp;
424 for(long row=0; row<r; ++row)
425 {
426 for(long col=row; col<c; ++col)
427 {
428 if(row == col)
429 continue;
430 temp = getDirect(row,col);
431 setDirect(row,col, getDirect(col,row));
432 setDirect(col,row, temp);
433 }
434 }
435 }
436
437 void transpose(matrix& B)
438 {
439 if(!elements)
440 return;
441 B.create(c,r);
442 for(int _r=0; _r<r; ++_r)
443 {
444 for(int _c=0; _c<c; ++_c)
445 {
446 Complex c = (*this)(_r, _c);
447 B(_c, _r) = creal(c) - cimag(c)*I;
448 }
449 //B.set(_c,_r, get(_r,_c));
450 }
451 }
452
453 void subtract(Complex& c)
454 {
455 if(!elements)
456 return;
457
458 for(int i=0; i<n; ++i)
459 {
460 elements[i] = elements[i] = c;
461 }
462 }
463
464 void subtract(matrix* B, matrix* output)
465 {
466 if(!B || !output)
467 return;
468
469 if(!elements || !B->getElements() || getRows() != B->getRows() || getCols() != B->getCols())
470 {
471 throw std::string("Error - cannot subtract matrix\n");
472 return;
473 }
474
475 output->create(getRows(), getCols());
476
477 for(int i=0; i<n; ++i)
478 {
479 output->elements[i] = elements[i] - B->getElements()[i];
480 }
481 }
482
483 void add(matrix* B, matrix* output)
484 {
485 if(!B || !output)
486 return;
487
488 if(!elements || !B->getElements() || getRows() != B->getRows() || getCols() != B->getCols())
489 {
490 throw std::string("Error - cannot add matrix\n");
491 return;
492 }
493
494 output->create(getRows(), getCols());
495
496 for(int i=0; i<n; ++i)
497 {
498 output->elements[i] = elements[i] + B->getElements()[i];
499 }
500 }
501
502 void random(double smallValue, double largeValue);
503 void add(matrix* B)
504 {
505 matrix* output = this;
506 if(!B)
507 return;
508
509 if(!elements || !B->getElements() || getRows() != B->getRows() || getCols() != B->getCols())
510 {
511 throw std::string("Error - cannot add matrix\n");
512 return;
513 }
514
515 //output->create(getRows(), getCols());
516
517 for(int i=0; i<n; ++i)
518 {
519 output->elements[i] = elements[i] + B->getElements()[i];
520 }
521 }
522
523 void dot(matrix* B, Complex& output); // returns this <dot> B
524 // set Trans* = 'T' for *^T
525 bool multiply(matrix& B, matrix& output, char TRANSB='N');
526
527 void tensor(matrix* B, matrix* output);
528
529 int getRows()
530 {
531 return r;
532 }
533 int getCols()
534 {
535 return c;
536 }
537 int getSize()
538 {
539 return n;
540 }
541 private:
542 int alg_min(int a, int b)
543 {
544 return (a < b) ? a : b;
545 }
546 int alg_max(int a, int b)
547 {
548 return (a > b) ? a : b;
549 }
550 Complex* elements;
551 int r, c, n;
552
553 bool noCleanup;
554
555 std::string varName;
556 };
557
558 typedef matrix mat;
559};
560
561#endif