Skip to content
Prev 999 / 10988 Next

[Rcpp-devel] Multiplication of ComplexVector?

Thank you for your interest in Rcpp.

Le 17/08/10 02:19, Christian Gunning a ?crit :
All is needed are binary operators on Rcomplex, I will add them in Rcpp 
before the next release. Or maybe someone else will follow my footsteps 
and do it.

One other thing about your code is that your loop is wrong :

for (j = 0; j++; j<n) {

should be :

for (j = 0; j<n; j++) {



Anyway, with some help, it works :

require( inline )
require( Rcpp )

inc <- '
Rcomplex operator*( const Rcomplex& lhs, const Rcomplex& rhs){
	Rcomplex y ;
	y.r = lhs.r * rhs.r - lhs.i * rhs.i ;
	y.i = lhs.r * rhs.i + rhs.r * lhs.i ;
	return y ;
}
Rcomplex operator+( const Rcomplex& lhs, const Rcomplex& rhs){
	Rcomplex y ;
	y.r = lhs.r + rhs.r ;
	y.i = lhs.i + rhs.i ;
	return y ;
}
// ...
'
fx <- cxxfunction( signature( i = "complex", ii = "complex" ), '
using namespace Rcpp ;
     ComplexVector y1(i), y2(ii);
     int n = y1.size() ;
     ComplexVector y3(n), y4(n);
     y3 = y1 * y2;  // no operator
     for (int j = 0; j<n ; j++) {
         y4[j] = y1[j] + y2[j]; // no operator
     }
     List z            = List::create(
     	_["*"] = y3,
     	_["+"] = y4,
     	_["y1*y2 + y1"] = y1 * y2 + y1,
     	_["y1 + y2*y2"] = y1 + y2 * y2
     	) ;
     return z ;
', plugin = "Rcpp", includes = inc )

y1 <- 1:10*(1+1i)
y2 <- 2*(1:10)*(1+1i)
stopifnot( identical(
	fx( y1, y2 ),
	list( "*" = y1 * y2, "+" = y1 + y2,
		"y1*y2 + y1" = y1 * y2 + y1,
		"y1 + y2*y2" = y1 + y2 * y2
	) ) )



 > fx( y1, y2 )
$`*`
  [1] 0+  4i 0+ 16i 0+ 36i 0+ 64i 0+100i 0+144i 0+196i 0+256i 0+324i 0+400i

$`+`
  [1]  3+ 3i  6+ 6i  9+ 9i 12+12i 15+15i 18+18i 21+21i 24+24i 27+27i 30+30i

$`y1*y2 + y1`
  [1]  1+  5i  2+ 18i  3+ 39i  4+ 68i  5+105i  6+150i  7+203i  8+264i 
9+333i
[10] 10+410i

$`y1 + y2*y2`
  [1]  1+  9i  2+ 34i  3+ 75i  4+132i  5+205i  6+294i  7+399i  8+520i 
9+657i
[10] 10+810i


Romain

  
    

Thread (24 messages)