b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2util.H
Go to the documentation of this file.
1//------------------------------------------------------------------------
2// b2util.H --
3//
4//
5// written by Mathias Doreille
6//
7// Copyright (c) 2004-2012,2017
8// SMR Engineering & Development SA
9// 2502 Bienne, Switzerland
10//
11// All Rights Reserved. Proprietary source code. The contents of
12// this file may not be disclosed to third parties, copied or
13// duplicated in any form, in whole or in part, without the prior
14// written permission of SMR.
15//------------------------------------------------------------------------
16
17#ifndef __B2UTIL_H__
18#define __B2UTIL_H__
19
20#include <algorithm>
21#include <array>
22#include <cerrno>
23#include <complex>
24#include <cstdlib>
25#include <cstring>
26#include <functional>
27#include <ostream>
28#include <set>
29#include <string>
30#include <typeinfo>
31#include <vector>
32
33#include "b2ppconfig.h"
34
35#if HAVE_GCC_ABI_DEMANGLE
36#include <cxxabi.h>
37#endif
38
39#include "b2csda.H"
40
46namespace b2000 {
47
48using std::size_t;
49
54class StringList : public std::vector<std::string> {
55public:
56 StringList() {}
57 StringList(const std::string& a1) { push_back(a1); }
58 StringList(const std::string& a1, const std::string& a2) {
59 push_back(a1);
60 push_back(a2);
61 }
62 StringList(const std::string& a1, const std::string& a2, const std::string& a3) {
63 push_back(a1);
64 push_back(a2);
65 push_back(a3);
66 }
67 StringList(
68 const std::string& a1, const std::string& a2, const std::string& a3,
69 const std::string& a4) {
70 push_back(a1);
71 push_back(a2);
72 push_back(a3);
73 push_back(a4);
74 }
75 StringList(
76 const std::string& a1, const std::string& a2, const std::string& a3,
77 const std::string& a4, const std::string& a5) {
78 push_back(a1);
79 push_back(a2);
80 push_back(a3);
81 push_back(a4);
82 push_back(a5);
83 }
84 StringList(
85 const std::string& a1, const std::string& a2, const std::string& a3,
86 const std::string& a4, const std::string& a5, const std::string& a6) {
87 push_back(a1);
88 push_back(a2);
89 push_back(a3);
90 push_back(a4);
91 push_back(a5);
92 push_back(a6);
93 }
94 StringList(
95 const std::string& a1, const std::string& a2, const std::string& a3,
96 const std::string& a4, const std::string& a5, const std::string& a6,
97 const std::string& a7) {
98 push_back(a1);
99 push_back(a2);
100 push_back(a3);
101 push_back(a4);
102 push_back(a5);
103 push_back(a6);
104 push_back(a7);
105 }
106 StringList(
107 const std::string& a1, const std::string& a2, const std::string& a3,
108 const std::string& a4, const std::string& a5, const std::string& a6,
109 const std::string& a7, const std::string& a8) {
110 push_back(a1);
111 push_back(a2);
112 push_back(a3);
113 push_back(a4);
114 push_back(a5);
115 push_back(a6);
116 push_back(a7);
117 push_back(a8);
118 }
119 StringList(
120 const std::string& a1, const std::string& a2, const std::string& a3,
121 const std::string& a4, const std::string& a5, const std::string& a6,
122 const std::string& a7, const std::string& a8, const std::string& a9) {
123 push_back(a1);
124 push_back(a2);
125 push_back(a3);
126 push_back(a4);
127 push_back(a5);
128 push_back(a6);
129 push_back(a7);
130 push_back(a8);
131 push_back(a9);
132 }
133 StringList(
134 const std::string& a1, const std::string& a2, const std::string& a3,
135 const std::string& a4, const std::string& a5, const std::string& a6,
136 const std::string& a7, const std::string& a8, const std::string& a9,
137 const std::string& a10) {
138 push_back(a1);
139 push_back(a2);
140 push_back(a3);
141 push_back(a4);
142 push_back(a5);
143 push_back(a6);
144 push_back(a7);
145 push_back(a8);
146 push_back(a9);
147 push_back(a10);
148 }
149 StringList(
150 const std::string& a1, const std::string& a2, const std::string& a3,
151 const std::string& a4, const std::string& a5, const std::string& a6,
152 const std::string& a7, const std::string& a8, const std::string& a9,
153 const std::string& a10, const std::string& a11) {
154 push_back(a1);
155 push_back(a2);
156 push_back(a3);
157 push_back(a4);
158 push_back(a5);
159 push_back(a6);
160 push_back(a7);
161 push_back(a8);
162 push_back(a9);
163 push_back(a10);
164 push_back(a11);
165 }
166 StringList(
167 const std::string& a1, const std::string& a2, const std::string& a3,
168 const std::string& a4, const std::string& a5, const std::string& a6,
169 const std::string& a7, const std::string& a8, const std::string& a9,
170 const std::string& a10, const std::string& a11, const std::string& a12) {
171 push_back(a1);
172 push_back(a2);
173 push_back(a3);
174 push_back(a4);
175 push_back(a5);
176 push_back(a6);
177 push_back(a7);
178 push_back(a8);
179 push_back(a9);
180 push_back(a10);
181 push_back(a11);
182 push_back(a12);
183 }
184 StringList(
185 const std::string& a1, const std::string& a2, const std::string& a3,
186 const std::string& a4, const std::string& a5, const std::string& a6,
187 const std::string& a7, const std::string& a8, const std::string& a9,
188 const std::string& a10, const std::string& a11, const std::string& a12,
189 const std::string& a13) {
190 push_back(a1);
191 push_back(a2);
192 push_back(a3);
193 push_back(a4);
194 push_back(a5);
195 push_back(a6);
196 push_back(a7);
197 push_back(a8);
198 push_back(a9);
199 push_back(a10);
200 push_back(a11);
201 push_back(a12);
202 push_back(a13);
203 }
204 StringList(
205 const std::string& a1, const std::string& a2, const std::string& a3,
206 const std::string& a4, const std::string& a5, const std::string& a6,
207 const std::string& a7, const std::string& a8, const std::string& a9,
208 const std::string& a10, const std::string& a11, const std::string& a12,
209 const std::string& a13, const std::string& a14) {
210 push_back(a1);
211 push_back(a2);
212 push_back(a3);
213 push_back(a4);
214 push_back(a5);
215 push_back(a6);
216 push_back(a7);
217 push_back(a8);
218 push_back(a9);
219 push_back(a10);
220 push_back(a11);
221 push_back(a12);
222 push_back(a13);
223 push_back(a14);
224 }
225 StringList(
226 const std::string& a1, const std::string& a2, const std::string& a3,
227 const std::string& a4, const std::string& a5, const std::string& a6,
228 const std::string& a7, const std::string& a8, const std::string& a9,
229 const std::string& a10, const std::string& a11, const std::string& a12,
230 const std::string& a13, const std::string& a14, const std::string& a15) {
231 push_back(a1);
232 push_back(a2);
233 push_back(a3);
234 push_back(a4);
235 push_back(a5);
236 push_back(a6);
237 push_back(a7);
238 push_back(a8);
239 push_back(a9);
240 push_back(a10);
241 push_back(a11);
242 push_back(a12);
243 push_back(a13);
244 push_back(a14);
245 push_back(a15);
246 }
247 StringList(
248 const std::string& a1, const std::string& a2, const std::string& a3,
249 const std::string& a4, const std::string& a5, const std::string& a6,
250 const std::string& a7, const std::string& a8, const std::string& a9,
251 const std::string& a10, const std::string& a11, const std::string& a12,
252 const std::string& a13, const std::string& a14, const std::string& a15,
253 const std::string& a16) {
254 push_back(a1);
255 push_back(a2);
256 push_back(a3);
257 push_back(a4);
258 push_back(a5);
259 push_back(a6);
260 push_back(a7);
261 push_back(a8);
262 push_back(a9);
263 push_back(a10);
264 push_back(a11);
265 push_back(a12);
266 push_back(a13);
267 push_back(a14);
268 push_back(a15);
269 push_back(a16);
270 }
271 StringList(
272 const std::string& a1, const std::string& a2, const std::string& a3,
273 const std::string& a4, const std::string& a5, const std::string& a6,
274 const std::string& a7, const std::string& a8, const std::string& a9,
275 const std::string& a10, const std::string& a11, const std::string& a12,
276 const std::string& a13, const std::string& a14, const std::string& a15,
277 const std::string& a16, const std::string& a17) {
278 push_back(a1);
279 push_back(a2);
280 push_back(a3);
281 push_back(a4);
282 push_back(a5);
283 push_back(a6);
284 push_back(a7);
285 push_back(a8);
286 push_back(a9);
287 push_back(a10);
288 push_back(a11);
289 push_back(a12);
290 push_back(a13);
291 push_back(a14);
292 push_back(a15);
293 push_back(a16);
294 push_back(a17);
295 }
296 StringList(
297 const std::string& a1, const std::string& a2, const std::string& a3,
298 const std::string& a4, const std::string& a5, const std::string& a6,
299 const std::string& a7, const std::string& a8, const std::string& a9,
300 const std::string& a10, const std::string& a11, const std::string& a12,
301 const std::string& a13, const std::string& a14, const std::string& a15,
302 const std::string& a16, const std::string& a17, const std::string& a18) {
303 push_back(a1);
304 push_back(a2);
305 push_back(a3);
306 push_back(a4);
307 push_back(a5);
308 push_back(a6);
309 push_back(a7);
310 push_back(a8);
311 push_back(a9);
312 push_back(a10);
313 push_back(a11);
314 push_back(a12);
315 push_back(a13);
316 push_back(a14);
317 push_back(a15);
318 push_back(a16);
319 push_back(a17);
320 push_back(a18);
321 }
322 StringList(
323 const std::string& a1, const std::string& a2, const std::string& a3,
324 const std::string& a4, const std::string& a5, const std::string& a6,
325 const std::string& a7, const std::string& a8, const std::string& a9,
326 const std::string& a10, const std::string& a11, const std::string& a12,
327 const std::string& a13, const std::string& a14, const std::string& a15,
328 const std::string& a16, const std::string& a17, const std::string& a18,
329 const std::string& a19) {
330 push_back(a1);
331 push_back(a2);
332 push_back(a3);
333 push_back(a4);
334 push_back(a5);
335 push_back(a6);
336 push_back(a7);
337 push_back(a8);
338 push_back(a9);
339 push_back(a10);
340 push_back(a11);
341 push_back(a12);
342 push_back(a13);
343 push_back(a14);
344 push_back(a15);
345 push_back(a16);
346 push_back(a17);
347 push_back(a18);
348 push_back(a19);
349 }
350 StringList(
351 const std::string& a1, const std::string& a2, const std::string& a3,
352 const std::string& a4, const std::string& a5, const std::string& a6,
353 const std::string& a7, const std::string& a8, const std::string& a9,
354 const std::string& a10, const std::string& a11, const std::string& a12,
355 const std::string& a13, const std::string& a14, const std::string& a15,
356 const std::string& a16, const std::string& a17, const std::string& a18,
357 const std::string& a19, const std::string& a20) {
358 push_back(a1);
359 push_back(a2);
360 push_back(a3);
361 push_back(a4);
362 push_back(a5);
363 push_back(a6);
364 push_back(a7);
365 push_back(a8);
366 push_back(a9);
367 push_back(a10);
368 push_back(a11);
369 push_back(a12);
370 push_back(a13);
371 push_back(a14);
372 push_back(a15);
373 push_back(a16);
374 push_back(a17);
375 push_back(a18);
376 push_back(a19);
377 push_back(a20);
378 }
379};
380
381inline std::string demangle(const std::type_info& t) {
382#if HAVE_GCC_ABI_DEMANGLE
383 int status = 0;
384 char* c_name = abi::__cxa_demangle(t.name(), 0, 0, &status);
385 std::string res(c_name);
386 ::free(c_name);
387 return res;
388#else
389 return t.name();
390#endif
391}
392
393inline std::ostream& operator<<(std::ostream& out, const std::type_info& t) {
394#if HAVE_GCC_ABI_DEMANGLE
395 int status = 0;
396 char* c_name = abi::__cxa_demangle(t.name(), 0, 0, &status);
397 out << c_name;
398 ::free(c_name);
399 return out;
400#else
401 out << t.name();
402 return out;
403#endif
404}
405
406inline std::string strerror(const int error) {
407 std::array<char, 256> buf{};
408
409#ifdef _WIN32
410 if (!strerror_s(buf.data(), buf.size(), error)) { return buf.data(); }
411#elif (_POSIX_C_SOURCE >= 200112L) && !_GNU_SOURCE && !__ve
412 if(!strerror_r(error, buf.data(), buf.size())
413 return buf.data();
414#elif __NEC__ && __ve
415 // The NEC vector compiler does not know strerror_r, just use
416 // a the error message without trying to look up the error code.
417 return std::string{"unknown error " + std::to_string(error)};
418#else
419 const char* msg = strerror_r(error, buf.data(), buf.size());
420 if (msg != nullptr) { return msg; }
421#endif
422
423 return std::string{"unknown error " + std::to_string(error)};
424}
425
426struct CompareFirstOfPair {
427 template <typename T1, typename T2>
428 bool operator()(const T1& a, const T2& b) const {
429 return a.first < b.first;
430 }
431};
432
433template <typename T1, typename T2, typename T3, typename T4>
434std::pair<T1, T2>& operator+=(std::pair<T1, T2>& i, const std::pair<T3, T4>& j) {
435 i.first += j.first;
436 i.second += j.second;
437 return i;
438}
439
440inline double norm(const float& v) { return std::abs(v); }
441
442inline double norm(const double& v) { return std::abs(v); }
443
444template <typename T>
445inline b2000::csda<T> norm(const b2000::csda<T>& v) {
446 return std::abs(v);
447}
448
449template <typename T>
450inline T norm(const std::complex<T>& v) {
451 return std::abs(v.real()) + std::abs(v.imag());
452}
453
454template <typename T>
455inline T max_abs(const T& a, const T& b) {
456 return std::max(b2000::norm(a), b2000::norm(b));
457}
458
459template <typename T>
460inline b2000::csda<T> max_abs(const b2000::csda<T>& a, const b2000::csda<T>& b) {
461 return std::max(b2000::norm(a), b2000::norm(b));
462}
463
464template <typename T>
465inline std::complex<T> max_abs(const std::complex<T>& a, const std::complex<T>& b) {
466 const T m = std::max(b2000::norm(a), b2000::norm(b));
467 return std::complex<T>(m);
468}
469
470inline double real(double v) { return v; }
471
472inline double imag(double v) { return 0; }
473
474template <typename T>
475inline T real(b2000::csda<T> v) {
476 return v.real();
477}
478
479template <typename T>
480inline T imag(b2000::csda<T> v) {
481 return v.imag();
482}
483
484template <typename T>
485inline T real(std::complex<T> v) {
486 return v.real();
487}
488
489template <typename T>
490inline T imag(std::complex<T> v) {
491 return v.imag();
492}
493
494inline bool is_invalid_sqrt(const float v) { return v < 0; }
495
496inline bool is_invalid_sqrt(const double v) { return v < 0; }
497
498template <typename T>
499inline bool is_invalid_sqrt(const b2000::csda<T> v) {
500 return false;
501}
502
503template <typename T>
504inline bool is_invalid_sqrt(const std::complex<T> v) {
505 return false;
506}
507
508template <typename _Tp1>
509struct auto_ptr_array_ref {
510 _Tp1* _M_ptr;
511
512 explicit auto_ptr_array_ref(_Tp1* __p) : _M_ptr(__p) {}
513};
514
519template <typename _Tp>
520class auto_ptr_array {
521private:
522 _Tp* _M_ptr;
523
524public:
525 typedef _Tp element_type;
526
527 explicit auto_ptr_array(element_type* __p = 0) noexcept : _M_ptr(__p) {}
528
529 auto_ptr_array(auto_ptr_array& __a) noexcept : _M_ptr(__a.release()) {}
530
531 template <typename _Tp1>
532 auto_ptr_array(auto_ptr_array<_Tp1>& __a) noexcept : _M_ptr(__a.release()) {}
533
534 auto_ptr_array& operator=(auto_ptr_array& __a) noexcept {
535 reset(__a.release());
536 return *this;
537 }
538
539 template <typename _Tp1>
540 auto_ptr_array& operator=(auto_ptr_array<_Tp1>& __a) noexcept {
541 reset(__a.release());
542 return *this;
543 }
544
545 ~auto_ptr_array() { delete[] _M_ptr; }
546
547 element_type& operator*() const noexcept { return *_M_ptr; }
548
549 element_type* operator->() const noexcept { return _M_ptr; }
550
551 element_type& operator[](size_t i) const noexcept { return _M_ptr[i]; }
552
553 element_type* get() const noexcept { return _M_ptr; }
554
555 element_type* release() noexcept {
556 element_type* __tmp = _M_ptr;
557 _M_ptr = 0;
558 return __tmp;
559 }
560
561 void reset(element_type* __p = 0) noexcept {
562 if (__p != _M_ptr) {
563 delete[] _M_ptr;
564 _M_ptr = __p;
565 }
566 }
567
568 auto_ptr_array(auto_ptr_array_ref<element_type> __ref) noexcept : _M_ptr(__ref._M_ptr) {}
569
570 auto_ptr_array& operator=(auto_ptr_array_ref<element_type> __ref) noexcept {
571 if (__ref._M_ptr != this->get()) {
572 delete _M_ptr;
573 _M_ptr = __ref._M_ptr;
574 }
575 return *this;
576 }
577
578 template <typename _Tp1>
579 operator auto_ptr_array_ref<_Tp1>() noexcept {
580 return auto_ptr_array_ref<_Tp1>(this->release());
581 }
582
583 template <typename _Tp1>
584 operator auto_ptr_array<_Tp1>() noexcept {
585 return auto_ptr_array<_Tp1>(this->release());
586 }
587};
588
589template <class T1, class T2>
590struct pair_iterator_reference {
591 T1& first;
592 T2& second;
593
594 pair_iterator_reference() : first(), second() {}
595
596 pair_iterator_reference(T1& a, T2& b) : first(a), second(b) {}
597
598 operator std::pair<T1, T2>() const { return std::pair<T1, T2>(first, second); }
599
600 pair_iterator_reference operator=(const pair_iterator_reference& p) {
601 first = p.first;
602 second = p.second;
603 return *this;
604 }
605
606 pair_iterator_reference operator=(const std::pair<T1, T2>& p) {
607 first = p.first;
608 second = p.second;
609 return *this;
610 }
611 friend void swap(pair_iterator_reference a, pair_iterator_reference b) {
612 using std::swap; // bring in swap for built-in types
613 swap(a.first, b.first);
614 swap(a.second, b.second);
615 }
616};
617
618template <typename Iterator1, typename Iterator2>
619struct pair_iterator {
620 using iterator_category = typename Iterator1::iterator_category;
621 using value_type = std::pair<typename Iterator1::value_type, typename Iterator2::value_type>;
622 using difference_type = typename Iterator1::difference_type;
623 using pointer = std::pair<typename Iterator1::pointer, typename Iterator2::pointer>;
624 using reference = pair_iterator_reference<typename Iterator1::value_type, typename Iterator2::value_type>;
625
626 Iterator1 first;
627 Iterator2 second;
628 pair_iterator() : first(), second() {}
629 pair_iterator(const Iterator1& a, const Iterator2& b) : first(a), second(b) {}
630 pair_iterator(const pair_iterator& i) : first(i.first), second(i.second) {}
631 pair_iterator& operator=(const pair_iterator& i) {
632 first = i.first;
633 second = i.second;
634 return *this;
635 }
636 reference operator*() { return reference(*first, *second); }
637 value_type operator*() const { return value_type(*first, *second); }
638 reference operator->() { return reference(*first, *second); }
639 pair_iterator& operator++() {
640 ++first;
641 ++second;
642 return *this;
643 }
644 pair_iterator& operator--() {
645 --first;
646 --second;
647 return *this;
648 }
649 pair_iterator operator++(int) {
650 pair_iterator tmp(*this);
651 ++first;
652 ++second;
653 return tmp;
654 }
655 pair_iterator operator--(int) {
656 pair_iterator tmp(*this);
657 --first;
658 --second;
659 return tmp;
660 }
661 pair_iterator& operator+=(const difference_type n) {
662 first += n;
663 second += n;
664 return *this;
665 }
666 pair_iterator& operator-=(const difference_type n) {
667 first -= n;
668 second -= n;
669 return *this;
670 }
671 pair_iterator operator+(const difference_type n) const {
672 return pair_iterator(first + n, second + n);
673 }
674 friend pair_iterator operator+(const difference_type n, const pair_iterator i) {
675 return pair_iterator(i.first + n, i.second + n);
676 }
677 pair_iterator operator-(const difference_type n) const {
678 return pair_iterator(first - n, second - n);
679 }
680 difference_type operator-(const pair_iterator i) const { return first - i.first; }
681 reference operator[](const difference_type n) { return reference(first[n], second[n]); }
682 bool operator==(const pair_iterator i) const { return first == i.first; }
683 bool operator!=(const pair_iterator i) const { return first != i.first; }
684 bool operator<(const pair_iterator i) const { return first < i.first; }
685 bool operator>(const pair_iterator i) const { return first > i.first; }
686 bool operator<=(const pair_iterator i) const { return first <= i.first; }
687 bool operator>=(const pair_iterator i) const { return first >= i.first; }
688};
689
690template <class Key, class Compare = std::less<Key> >
691class vector_set : public std::vector<Key> {
692public:
693 vector_set(const Compare& comp = Compare()) : comp_o(comp), set_o(Comp(*this, comp_o)) {
694 current = set_o.end();
695 }
696
697 size_t insert(const Key& key) {
698 if (current != set_o.end()) {
699 if (comp_o.equal(key, (*this)[*current])) {
700 return *current;
701 } else {
702 ++current;
703 }
704 }
705 const size_t pos = std::vector<Key>::size();
706 std::vector<Key>::push_back(key);
707 current = set_o.insert(current, pos);
708 if (*current != pos) { std::vector<Key>::pop_back(); }
709 return *current;
710 }
711
712 size_t insert_cycle(const Key& key) {
713 if (current != set_o.end()) {
714 if (comp_o.equal(key, (*this)[*current])) { return *current; }
715 ++current;
716 if (current == set_o.end()) { current = set_o.begin(); }
717 if (comp_o.equal(key, (*this)[*current])) { return *current; }
718 }
719 const size_t pos = std::vector<Key>::size();
720 std::vector<Key>::push_back(key);
721 current = set_o.insert(current, pos);
722 if (*current != pos) { std::vector<Key>::pop_back(); }
723 return *current;
724 }
725
726 size_t insert_equal_or_upper_to_top(const Key& key) {
727 if (!std::vector<Key>::empty() && comp_o.equal(key, std::vector<Key>::back())) {
728 return std::vector<Key>::size() - 1;
729 }
730 std::vector<Key>::push_back(key);
731 return std::vector<Key>::size() - 1;
732 }
733
734private:
735 Compare comp_o;
736 struct Comp : public std::function<bool(size_t, size_t)> {
737 Comp(const std::vector<Key>& vect_, const Compare& comp_o_)
738 : vect(vect_), comp_o(comp_o_) {}
739 const std::vector<Key>& vect;
740 Compare comp_o;
741 constexpr bool operator()(const size_t& a, const size_t& b) const {
742 return comp_o(vect[a], vect[b]);
743 }
744 };
745 typedef std::set<size_t, Comp> set_t;
746 set_t set_o;
747 typename set_t::iterator current;
748};
749
769template <typename T>
770inline void solve_cubic_equation(
771 const T a, // coefficient of x^3
772 const T b, // coefficient of x^2
773 const T c, // coefficient of x
774 const T d, // constant term
775 int& nsol, // # of distinct solutions
776 T& x1, T& x2, T& x3) {
777 if (a == 0) {
778 if (b == 0) {
779 if (c == 0) {
780 if (d == 0) {
781 nsol = 1;
782 x1 = 0;
783 } else {
784 nsol = 0;
785 }
786 } else {
787 nsol = 1;
788 x1 = -d / c;
789 }
790 } else {
791 T delta = c * c - 4 * b * d;
792 if (delta < 0) {
793 nsol = 0;
794 return;
795 }
796 delta = std::sqrt(delta);
797 x1 = (-c + delta) / (2 * b);
798 if (delta == 0) {
799 nsol = 1;
800 } else {
801 x2 = (-c - delta) / (2 * b);
802 nsol = 2;
803 }
804 }
805 } else {
806 const T a1 = b / a;
807 const T a2 = c / a;
808 const T a3 = d / a;
809 const T Q = (a1 * a1 - 3.0 * a2) / 9.0;
810 const T R = (2.0 * a1 * a1 * a1 - 9.0 * a1 * a2 + 27.0 * a3) / 54.0;
811 const T R2_Q3 = R * R - Q * Q * Q;
812
813 if (R2_Q3 <= 0) {
814 nsol = 3;
815 const T theta = std::acos(R / std::sqrt(Q * Q * Q));
816 const T pi = std::acos(T(-1));
817 x1 = -2.0 * std::sqrt(Q) * std::cos(theta / 3.0) - a1 / 3.0;
818 x2 = -2.0 * std::sqrt(Q) * std::cos((theta + 2.0 * pi) / 3.0) - a1 / 3.0;
819 x3 = -2.0 * std::sqrt(Q) * std::cos((theta + 4.0 * pi) / 3.0) - a1 / 3.0;
820 } else {
821 nsol = 1;
822 x1 = std::pow(std::sqrt(R2_Q3) + std::abs(R), 1 / 3.0);
823 x1 += Q / x1;
824 x1 *= (R < 0.0) ? 1.0 : -1.0;
825 x1 -= a1 / 3.0;
826 }
827 }
828}
829
830// small limited shared_ptr implementation in the b2000 namespace, until we switch to c++-03
831template <typename T>
832class shared_ptr {
833public:
834 shared_ptr() : ptr_(0), ref_count_(nullptr) {}
835
836 shared_ptr(T* p) : ptr_(p), ref_count_(p ? new int : 0) { inc_ref(); }
837
838 shared_ptr(const shared_ptr& rhs) : ptr_(rhs.ptr_), ref_count_(rhs.ref_count_) { inc_ref(); }
839
840 ~shared_ptr() {
841 if (ref_count_ && dec_ref() == 0) {
842 delete ptr_;
843 delete ref_count_;
844 }
845 }
846
847 T* get() { return ptr_; }
848
849 const T* get() const { return ptr_; }
850
851 void swap(shared_ptr& rhs) {
852 std::swap(ptr_, rhs.ptr_);
853 std::swap(ref_count_, rhs.ref_count_);
854 }
855
856 shared_ptr& operator=(const shared_ptr& rhs) {
857 shared_ptr tmp(rhs);
858 this->swap(tmp);
859 return *this;
860 }
861
862 T* operator->() const { return ptr_; }
863
864 T* operator*() const { return ptr_; }
865
866private:
867 void inc_ref() {
868 if (ref_count_) { ++(*ref_count_); }
869 }
870
871 int dec_ref() { return --(*ref_count_); }
872
873 T* ptr_;
874 int* ref_count_;
875};
876
877template <typename T, int i1>
878inline void fill_array(T (&t)[i1], const T& v) {
879 std::fill_n(&t[0], i1, v);
880}
881
882template <typename T, int i1, int i2>
883inline void fill_array(T (&t)[i1][i2], const T& v) {
884 std::fill_n(&t[0][0], i1 * i2, v);
885}
886
887template <typename T, int i1, int i2, int i3>
888inline void fill_array(T (&t)[i1][i2][i3], const T& v) {
889 std::fill_n(&t[0][0][0], i1 * i2 * i3, v);
890}
891
892template <typename T, int i1, int i2, int i3, int i4>
893inline void fill_array(T (&t)[i1][i2][i3][i4], const T& v) {
894 std::fill_n(&t[0][0][0][0], i1 * i2 * i3 * i4, v);
895}
896
897template <typename T, int i1>
898inline void copy_array(const T (&t)[i1], T (&d)[i1]) {
899 std::copy(&t[0], &t[0] + i1, &d[0]);
900}
901
902template <typename T, int i1, int i2>
903inline void copy_array(const T (&t)[i1][i2], T (&d)[i1][i2]) {
904 std::copy(&t[0][0], &t[0][0] + i1 * i2, &d[0][0]);
905}
906
907template <typename T, int i1, int i2, int i3>
908inline void copy_array(const T (&t)[i1][i2][i3], T (&d)[i1][i2][i3]) {
909 std::copy(&t[0][0][0], &t[0][0][0] + i1 * i2 * i3, &d[0][0][0]);
910}
911
912template <typename T, int i1, int i2, int i3, int i4>
913inline void copy_array(const T (&t)[i1][i2][i3][i4], T (&d)[i1][i2][i3][i4]) {
914 std::copy(&t[0][0][0][0], &t[0][0][0][0] + i1 * i2 * i3 * i4, &d[0][0][0][0]);
915}
916} // namespace b2000
917
918#endif
b2000::details::operator>=
bool operator>=(const csda< T1 > &a, const csda< T2 > &b)
Comparison of two csda numbers is performed on the real part only.
Definition b2csda.H:316
b2000::details::operator<=
bool operator<=(const csda< T1 > &a, const csda< T2 > &b)
Comparison of two csda numbers is performed on the real part only.
Definition b2csda.H:298
b2000::details::norm
csda< T > norm(const csda< T > &a)
Definition b2csda.H:343
b2000::details::operator==
bool operator==(const csda< T1 > &a, const csda< T2 > &b)
Comparison of two csda numbers is performed on the real part only.
Definition b2csda.H:226
b2000::details::operator>
bool operator>(const csda< T1 > &a, const csda< T2 > &b)
Comparison of two csda numbers is performed on the real part only.
Definition b2csda.H:280
b2000::details::operator!=
bool operator!=(const csda< T1 > &a, const csda< T2 > &b)
Comparison of two csda numbers is performed on the real part only.
Definition b2csda.H:244
b2000::details::operator<
bool operator<(const csda< T1 > &a, const csda< T2 > &b)
Comparison of two csda numbers is performed on the real part only.
Definition b2csda.H:262
b2000::StringList
Definition b2util.H:54
b2000::auto_ptr_array
Definition b2util.H:520
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::solve_cubic_equation
void solve_cubic_equation(const T a, const T b, const T c, const T d, int &nsol, T &x1, T &x2, T &x3)
Definition b2util.H:770