I did not really expect the following example to work, but indeed it does (g++ 4.6.4, with --std=c++0x):
#include <boost/multiprecision/float128.hpp>
#include <blitz/array.h>
#include <fftw3.h>
int main(int /*argc*/, char** /*argv*/)
{
//these are the same
std::cout << sizeof(std::complex<boost::multiprecision::float128>) << " " << sizeof(fftwq_complex) << std::endl;
typedef std::vector< std::complex<boost::multiprecision::float128> > boost128cvec;
//typedef std::vector<std::complex<boost::multiprecision::float128> , fftw::allocator< std::complex<boost::multiprecision::float128> > > boost128cvec;
//declare a std::vector consisting of std::complex<boost::multiprecision::float128>
boost128cvec test_vector3(12);
//casting its data storatge to fftwq_complex*
fftwq_complex* test_ptr3 = reinterpret_cast<fftwq_complex*>(test_vector3.data());
//also create a view to the same data as a blitz::Array
blitz::Array<std::complex<boost::multiprecision::float128>, 1> test_array3(test_vector3.data(), blitz::TinyVector<int, 1>(12), blitz::neverDeleteData);
test_vector3[3] = std::complex<boost::multiprecision::float128>(1.23,4.56);
//this line would not work with std::vector
test_array3 = sin(test_array3);
//this line would not work with the built-in type __float128
test_vector3[4] = sin(test_vector3[3]);
//all of those print the same numbers
std::cout << "fftw::vector: " << test_vector3[3].real() << " + i " << test_vector3[3].imag() << std::endl;
std::cout << "fftw_complex: " << (long double)test_ptr3[3][0] << " + i " << (long double)test_ptr3[3][1] << std::endl;
std::cout << "blitz: " << test_array3(3).real() << " + i " << test_array3(3).imag() << std::endl << std::endl;
}
Two remarks:
fftw
and blitz::Array
operations on the same data without the need to copy them around while at the same time being able to use generic funcionst like sin()
also for complex variables with quad precisionblitz
-part works fine, which is expected. But the surprise (to me) was, that the fftwq_complex*
part also works fine.fftw::allocator
is a simple replacement to std::allocator
which will use fftwq_malloc
to assure correct simd alignment, but that is not important for this question, so I left it out (at least I think that this is not important for this question)My Question is: How thin is the ice I'm stepping on?
You're pretty much save:
std::vector
is compatible with a C array (you can access a pointer to the first element via vector.data()
, as answered in this questionstd::complex<T>
is designed to be compatible with a Array of form T[2]
, which is compatible with FFTW. This is described in the FFTW documentation
C++ has its own complex template class, defined in the standard header file. Reportedly, the C++ standards committee has recently agreed to mandate that the storage format used for this type be binary-compatible with the C99 type, i.e. an array T[2] with consecutive real [0] and imaginary [1] parts. (See report http://www.open-std.org/jtc1/sc22/WG21/docs/papers/2002/n1388.pdf WG21/N1388.) Although not part of the official standard as of this writing, the proposal stated that: “This solution has been tested with all current major implementations of the standard library and shown to be working.” To the extent that this is true, if you have a variable complex *x, you can pass it directly to FFTW via reinterpret_cast(x).
The only thing to keep in mind is that the data()
gets invalidated if you add values to your vector.
boost::multiprecision::float128
and __float128
. The boost documentation gives no guarantee about this.
What can be done however, is to add some static asserts in your code, which fails if the conversion is not possible. This could look like this:
static_assert(std::is_standard_layout<float128>::value,"no standard type");
static_assert(sizeof(float128) == sizeof(__float128),"size mismatch");
Where sizeof
guarantees the same size of the boost type and __float128, and is_standard_layout checks that:
A pointer to a standard-layout class may be converted (with reinterpret_cast) to a pointer to its first non-static data member and vice versa.
Of course, this only gives hints if it works in the end, as you cannot say if the type is really a __float128
, but ab boost states their type is a thin wrapper around it, it should be fine. If their are changes in design or structure of float128
, the static assertions should fail.