[Rcpp-devel] Multiplication of ComplexVector?

[Romain Francois]

Thank you for your interest in Rcpp.

Le 17/08/10 02:19, Christian Gunning a écrit :
> Dear list,
>
> I'm trying to use the ComplexVector class, and I'm having difficulty
> with multiplication. I sugar multiplication, element-by-element
> multiplication (no * operator for either), and using temporary
> Rcomplex variables (don't really understand what I'm doing here).

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

> I was able to successfully import get("*") as an Rcpp::Function and do
> the multiplication through R. Is this the best option for now?
>
> Thanks,
> Christian Gunning
> University of New Mexico
>
> // doesn't work
> SEXP rcpp_hello_world(SEXP i, SEXP ii){
>      using namespace Rcpp ;
>
>      ComplexVector y1(i), y2(ii);
>      int n(y1.size()), j;
>      ComplexVector y3(n), y4(n);
>      y3 = y1 * y2;  // no operator
>      for (j = 0; j++; j<n) {
>          y4[j] = y1[j] + y2[j]; // no operator
>      }
>      List z            = List::create( y1, y2, y3) ;
>      return z ;
> }
>
> // does work, doing multiplication in R
> SEXP mymult(SEXP input1, SEXP input2) {
>      Rcpp::Environment base("package:base");
>      Rcpp::Function funbase = base["get"];
>      SEXP funbase1 = funbase("*");
>      Rcpp::Function funmult(funbase1);
>      Rcpp::ComplexVector xxa = funmult(input1, input2);
>      return xxa;
> }
>
> SEXP rcpp_hello_world(SEXP i, SEXP ii){
>      using namespace Rcpp ;
>      ComplexVector y3( mymult(i, ii) );
>      return y3 ;
> }


-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://bit.ly/bzoWrs : Rcpp svn revision 2000
|- http://bit.ly/b8VNE2 : Rcpp at LondonR, oct 5th
`- http://bit.ly/aAyra4 : highlight 0.2-2