bug in spray multiplication

[stla]

Hello Robin,

There's a problem with the spray mutiplication here:

https://github.com/RobinHankin/spray/blob/master/src/spray_ops.cpp#L164C13-L164C33

After the Sout[vsum] += x1*x2, one should add Sout.erase(vsum) if Sout[vsum]==0.

[RobinHankin]

thanks for this! It's good to hear from you. As it stands, zero entries are discarded by prepare():

suppressMessages(library("spray"))
spray(diag(2),c(1,1))*spray(diag(2),c(1,-1))
#>          val
#>  2 0  =    1
#>  0 2  =   -1

^{Created on 2024-04-05 with reprex v2.1.0}

I believe that your fix would not do any harm [at least, implementing it does not break any of the package tests], but is not needed because prepare() strips out the zeros. It was a long time ago now, but I do remember considering adding an erase statement along the lines you suggest, but rejected it on efficiency grounds: it strips out zero values for a key, but that key might subsequently have to be reinstated. This would happen in the example above.

Comments?

[stla]

Ah ok. But maybe in general one would do less operations in this way as compared to cleaning the whole final spray. No?

I'm currently working on ratio of polynomials in R, e.g.:

> x <- qlone(1)
> y <- qlone(2)
> z <- qlone(3)
> (y^4 + y^3 + y^2 + y + 1) / (x + 1)  *  (y - 1) / (x - 1)
[x^(0, 5) - x^()] / [x^(2) - x^()] 

The ultimate goal is to deal with polynomials whose coefficients are ratios of polynomials.

[RobinHankin]

Ah, it's not clear to me whether the burden of repeatedly erasing and recreating a value would outweigh the difficulty of cleaning the whole final spray (in both multiplication and addition); I guess it's application-specific. The cost of zero-removal is $O(n\log n)$, comparable to addition but much smaller than multiplication. The spray package has the [compile time] option to use either vector or deque containers, which might make a difference too.

...

I spent quite a bit of time trying to set up a class of objects where the coefficients of variable x were polynomials in variables other than x, the idea being that the expressions would be things like [1+y]x^2 + [1+4y^2]x + [3]x^0. This isn't quite as general as your idea of coefficients being ratios of other polynomials. I got bogged down in a recursive nightmare and gave up...but the work wasn't wasted because it eventually led to the mvp package. But I do like the idea of working with ratios of polynomials as coefficients.

Best wishes

Robin

[stla]

Yeah my comment was stupid.

My motivation for the rational polynomials as the coefficients is the following fact about the Jack polynomials: the coefficients of these polynomials depend on a parameter and from the definition these coefficients are fractions of polynomials in this parameter. However by an amazing theorem they are polynomial actually.

I already did the implementation in Haskell. Haskell is great for such stuff. I love this language.

[stla]

Talking about efficiency, you can improve the code for the power of a spray:

      Spray Result = unitSpray;
      while(n) {
        if(n & 1) {
          Result *= S;
        }
        n >>= 1;
        S *= S;
      }
      return Result;

[RobinHankin]

that is super clever, and I like it. I like it a lot. I'll need to overload the * operator to do it properly though.

[stla]

I'll need to overload the * operator to do it properly though.

If you want, you can take my code of my qspray package. I have a templated class Qspray. Your spray objects correspond to Qspray<double>.

[stla]

Seems to work :-)

library(symbolicQspray)
> 
> f <- function(a1, a2, X1, X2, X3) {
+   (a1/(a2^2+1)) * X1^2*X2  +  (a2+1) * X3
+ }
> 
> a1 <- qlone(1)
> a2 <- qlone(2)
> X1 <- Qlone(1)
> X2 <- Qlone(2)
> X3 <- Qlone(3)
> 
> ( Qspray <- f(a1, a2, X1, X2, X3) )
{ [ a1 ] %//% [ a2^2 + 1 ] } * X1^2.X2  +  { [a2 + 1] } * X3 
> Qspray^2
{ [ a1^2 ] %//% [ a2^4 + 2*a2^2 + 1 ] } * X1^4.X2^2  +  { [ 2*a1.a2 + 2*a1 ] %//% [ a2^2 + 1 ] } * X1^2.X2.X3  +  { [a2^2 + 2*a2 + 1] } * X3^2 
> Qspray - Qspray
0 
> 
> ( qspray <- evalSymbolicQspray(Qspray, a = c(2, 3)) )
1/5*X1^2.X2 + 4*X3 
> ( ratioOfQsprays <- evalSymbolicQspray(Qspray, X = c(4, 3, 2)) )
[ 48*a1 + 2*a2^3 + 2*a2^2 + 2*a2 + 2 ]  %//%  [ a2^2 + 1 ] 
> evalSymbolicQspray(Qspray, a = c(2, 3), X = c(4, 3, 2))
Big Rational ('bigq') :
[1] 88/5
> evalQspray(qspray, c(4, 3, 2))
Big Rational ('bigq') :
[1] 88/5
> evalRatioOfQsprays(ratioOfQsprays, c(2, 3))
Big Rational ('bigq') :
[1] 88/5

[stla]

It's possible that the code for the power performs one unnecessary multiplication. Fix:

      Spray Result = unitSpray;
      int b = 1; p = 0;
      while(n) {
        if(n & 1) {
          Result *= S;
          p += b;
          if(p+1 == n) break;
        }
        n >>= 1;
        S *= S;
        b *= 2;
      }
      return Result;

[stla]

Sorry, wrong. Correct:

      Spray Result = unitSpray;
      int n0 = n, b = 1; p = 0;
      while(n) {
        if(n & 1) {
          Result *= S;
          p += b;
          if(p == n0) break;
        }
        n >>= 1;
        S *= S;
        b *= 2;
      }
      return Result;