error in ValueMolienSeries

[dimpase]

ValueMolienSeries(m,0) must be 1 for any m. However, the following shows a bug:

gap> g:=SymplecticGroup(6,3);;
gap> h:=Stabilizer(g,Z(3)*[1,0,0,0,0,0]);;
gap> t:=CharacterTable(h);;
gap> chi:=Irr(t)[7];;
gap> chi[1];
9
gap> m:=MolienSeries(t,chi);;
gap> ValueMolienSeries(m,0);
26/27

I didn't look at the code yet, but I guess it computes an approximate value and gets it wrong.

[olexandr-konovalov]

Thanks - I can reproduce this both in the master branch and GAP 4.7.8.

[hulpke]

On Oct 22, 2015, at 8:27 AM, Dmitrii Pasechnik notifications@github.com wrote: I didn't look at the code yet, but I guess it computes an approximate value and gets it wrong.

No, as far as I know it calculates exactly, but i the past the Molien series code was a good way to note the existence (not necessarily find where…) of bugs in a broad range of code.

[olexandr-konovalov]

I've tried

g:=SymplecticGroup(6,3);;
h:=Stabilizer(g,Z(3)*[1,0,0,0,0,0]);;
t:=CharacterTable(h);;
chi:=Irr(t)[7];;
chi[1];
m:=MolienSeries(t,chi);;
ValueMolienSeries(m,0);

in all 4 combinations of 32/64-bit builds with/without GMP, and 26/27 is returned in all 4 cases.

[olexandr-konovalov]

I can reproduce this with:

• GAP 4.6.5
• GAP 4.5.7
• GAP 4.4.12.

[fingolfin]

But contrast this:

gap> IsRationalFunction(m);
true
gap> Value(m, [IndeterminateOfUnivariateRationalFunction(m)], [0]);
1

[fingolfin]

Asking Mathematica for a Taylor series expansion gives $1+2z^{12}+5z^{18}+...$, and we get this in GAP:

gap> List([0..20], i->ValueMolienSeries(m, i));
[ 26/27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 5, 0, 0 ]

So, that looks about right, except of course for the constant term, which the above call to Value also correctly computes.

In other words, this gives me some hope that the Molien series is computed correctly, and that the bug is instead in ValueMolienSeries. Perhaps in CoefficientTaylorSeries?

[fingolfin]

I just used Mathematica to cross check the results CoefficientTaylorSeries provides, and they seem to be 100% correct. Yet they sum up to -12130560, while the size of the table is 12597120, leading to the incorrect result.

[dimpase]

"they sum up to" - well, it's not converging series at 0, so I don't quite get what you mean. And how come you get '-' sign?

On 4 November 2015 at 09:09, Max Horn notifications@github.com wrote:

I just used Mathematica to cross check the results CoefficientTaylorSeries provides, and they seem to be 100% correct. Yet they sum up to -12130560, while the size of the table is 12597120, leading to the incorrect result.

— Reply to this email directly or view it on GitHub https://github.com/gap-system/gap/issues/300#issuecomment-153647167.

[dimpase]

First few coefficients of the series are certainly OK; it can be checked directly from the character table: here is how to check for degree up to 36, and it agrees with the Taylor series

   onedim:=Filtered(Irr(t),x->x[1]=1);
   gap> List([1..36],x->MatScalarProducts(onedim,SymmetricParts(t,[chi],x))[1][1]);
   [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 13, 0, 
     0, 3, 0, 0, 33, 0, 0, 15, 0, 0, 87 ]

And here is the expansion into the Taylor series at 0 up to degree 300, using Sage. Looks OK; e.g. terms of degrees not divisible by 3 are all 0, as predicted by the theory here.

37017617_z^300 + 402747472_z^297 + 372723885_z^294 + 342894227_z^291 + 316865230_z^288 + 290972607_z^285 + 268474186_z^282 + 246064266_z^279 + 226678408_z^276 + 207340813_z^273 + 190691644_z^270 + 174057606_z^267 + 159807508_z^264 + 145545869_z^261 + 133392728_z^258 + 121206905_z^255 + 110881127_z^252 + 100505949_z^249 + 91767098_z^246 + 82966767_z^243 + 75601497_z^240 + 68165881_z^237 + 61985315_z^234 + 55728595_z^231 + 50565794_z^228 + 45323689_z^225 + 41031791_z^222 + 36659580_z^219 + 33109962_z^216 + 29480472_z^213 + 26560457_z^210 + 23562749_z^207 + 21174366_z^204 + 18711343_z^201 + 16769706_z^198 + 14757312_z^195 + 13188959_z^192 + 11554406_z^189 + 10296251_z^186 + 8976871_z^183 + 7975064_z^180 + 6917011_z^177 + 6125505_z^174 + 5283100_z^171 + 4663042_z^168 + 3997230_z^165 + 3515934_z^162 + 2993905_z^159 + 2623960_z^156 + 2218148_z^153 + 1936809_z^150 + 1624206_z^147 + 1412787_z^144 + 1174275_z^141 + 1017371_z^138 + 837383_z^135 + 722562_z^132 + 588201_z^129 + 505499_z^126 + 406444_z^123 + 347888_z^120 + 275824_z^117 + 235143_z^114 + 183477_z^111 + 155869_z^108 + 119376_z^105 + 101075_z^102 + 75781_z^99 + 64007_z^96 + 46776_z^93 + 39474_z^90 + 27981_z^87 + 23632_z^84 + 16146_z^81 + 13683_z^78 + 8938_z^75 + 7651_z^72 + 4709_z^69 + 4091_z^66 + 2346_z^63 + 2099_z^60 + 1087_z^57 + 1022_z^54 + 461_z^51 + 472_z^48 + 178_z^45 + 203_z^42 + 58_z^39 + 87_z^36 + 15_z^33 + 33_z^30 + 3_z^27 + 13_z^24 + 5_z^18 + 2*z^12 + 1

[fingolfin]

"they sum up to" - well, it's not converging series at 1, so I don't quite get what you mean. And how come you get '-' sign?

"They" in my remark referred to the result of CoefficientTaylorSeries. So to understand my remark, you need to look at the code. But roughly summarized, as far as I understand it: The Molien series is represented internally as a sum of rational functions p(z) / (1-z^r)^k, divided by the table size, and a "correction factor" of +1 or -1. So, to compute the 0th term of the taylor expansion, GAP computes the 0th term of the taylor expansion of these summands (which is what CoefficientTaylorSeries does, and it does it correctly), then sums them up (which in this case yields -12130560), then multiplies with the correction factor (which in this case is -1, because the degree is 9, i.e. odd), and finally divides by the table size -- which apparently is 12597120.

Finally, evaluating ValueMolienSeries up to degree 300, and subtracing that from the expansion you give, yields -400000000*z^300+1/27 ... that's weird (did you have a specific reason to pick 300?).

[dimpase]

It seems that the summands of the form p(z) / (1-z^r)^k correspond to maximal cyclic subgroups, up to conjugation. I didn't know that their Molien series have such a nice form. One needs to write this down, or to find a printed source, before understanding this bug...

300 was picked by accident. And I missed the 1st character while cutting/pasting, so it should be 437017617*z^300. And thus the difference is in fact 1/27, as expected. In case you need more coefs:

1408938228_z^348 + 1313744182_z^345 + 1228179934_z^342 + 1143741754_z^339 + 1068049576_z^336 + 993307279_z^333 + 926495399_z^330 + 860477584_z^327 + 801638887_z^324 + 743456479_z^321 + 691761821_z^318 + 640603391_z^315 + 595296931_z^312 + 550421484_z^309 + 510815336_z^306 + 471548904_z^303 + 437017617*z^300 + ...

[fingolfin]

Yes, that fits with all the data: both the results of ValueMolienSeries, and the taylor expansion of the rational function m as computed by Mathematica. So, everything seems to be internally consistent, except for the coefficient at degree 0. And indeed, to debug this further, one needs a reference for the algorithm. The only hint for that I found is a reference in the manual entry for MolienSeries which points to "Neubüser, J., Pahlings, H. and Plesken, W. (Atkinson, M. D., Ed.), CAS; design and use of a system for the handling of characters of finite groups, in Computational group theory (Durham, 1982), Academic Press, London (1984), 195–247.", see http://www.ams.org/mathscinet-getitem?mr=760658 -- but I do not have access to that.

There is also http://www.gap-system.org/ForumArchive/Breuer.1/Thomas.1/Re__Moli.1/1.html

[dimpase]

I did a sanity check of CoefficientTaylorSeries by comparing its output (for the degree 0, i.e. for l==0) with the 1st coefficient of the numerator (numer), which should be equal, as we're evaluating coefficients of f(z)=numer(z)/(1-z^r)^k, and we just want f(0)=number(0). And there is no discrepancy. That is, the bug is not in the implementation of this function...

[olexandr-konovalov]

We have a copy of the Computational group theory (Durham, 1982) proceedings in St Andrews.

[fingolfin]

@dimpase Exactly what I trying to say :).

[dimpase]

I think I found the bug, sort of. Note that the rational function is formed OK in the reduced form (ratfun in MolienSeries), it's the computation of its expansion from stored summands is broken. And indeed, after the process of storing the summands and summing them up, in the loop starting at line 493 an extra term coming from something called pol gets added to the numerator of the rat.fun. in line 527

      numer:= numer + denom * UnivariatePolynomial( Rationals, pol );

and it happens to be a constant, equal to -Size(tbl)/27, which is exactly the discrepancy we see. But pol does not get added to the Molien series record, whereas it should be added to summands in some way, by right...

Needless to say, if pol had degree bigger than 0, it would affect coefficients of higher degree, so it's a miracle that this does not happen more often...

From reading the code, it appears that pol collects stuff from dividing p(z)'s by corresponding (1-z^r)^k's, so that the degrees of p(z)'s stay <= kr.

[dimpase]

I was looking at the old history, and saw

diff -r 9f5bd6c8ed70 -r 3bc3e5a127ad lib/ctblmoli.gi
--- a/lib/ctblmoli.gi   Tue May 10 08:30:45 2005 +0000
+++ b/lib/ctblmoli.gi   Tue Mar 09 15:19:52 2004 +0000
@@ -545,6 +545,10 @@
        #T avoid forming this quotient!
        SetIsUnivariateRationalFunction( series, true );

 +    # No polynomial is needed to correct the result
 +    # if everything worked well.
 +    Assert( 1, IsEmpty( pol ) );
 +
       # Set the info record.
       SetMolienSeriesInfo( series,
                      rec( summands    := numers,

that is to say, it was a sticking point before...

[olexandr-konovalov]

The commit message from May 10th, 2005 says removed a wrong assertion. The assertion has been added on July 7th, 1999 with commit message saying for documentation in ctblfuns.msk...

[dimpase]

right - the dates in my diff are backwards...

[olexandr-konovalov]

This example uses the 2-argument version of MolienSeries. Code coverage report shows that we test 1-argument and 3-argument versions, but never test 2-argument ones.

Just in case if this rings any bells, but I can not be sure.