Project Euler 486

[lavalamp]

So I've been doing these problems again lately and figured I'd take a look at the newest, number 486.

However, I think I must be missing something about the setup of it:

---------------
Let F₅(n) be the number of strings s such that:

s consists only '0's and '1's,
s has length at most n, and
s contains a palindromic substring of length at least 5.

For example, F₅(4) = 0, F₅(5) = 8, F₅(6) = 42 and F₅(11) = 3844.
---------------

Clearly F(4) = 0, because there are no substrings with length 5.

I'm also fairly happy with F(5) = 8.

My issue is with F(6) = 42. I can't figure how this is the case. For strings of length exactly 6, I count 28 with a required palindromic substring. Since the problem states that F(n) is defined as strings AT MOST length n, then we can include all of the length 5 strings too. So this would bring the answer for F(6) to 36.

Even if it were as simple as adding 1 extra bit to either the beginning or end of a palindromic string of 5, there would only be 4*8 = 32 solutions of length 6 (of course, some are double-ups that must be accounted for). But even in this case, F(6) = 40, still shy of 42.

I know that this is a new problem, so it could be a mistake. However 44 people have solved this one, so I'm unsure.

[Batalov]

Originally Posted by lavalamp

...s contains a palindromic substring of length at least 5. <-- (not just 5)
...

I had no problem with the formulation.
I reproduced D(5*10^9), too. I still have problem with final scale-up. I have an idea how to do it.

[lavalamp]

Ohhhhhhhhhhhh, I see. That makes it so much harder too.

[Batalov]

It is always interesting when one's answer to a problem is wrong not because the solution is wrong, but because of a bug in Pari.

A certain (rather straightforward) solution for PE 501 gives the wrong answer with Pari 2.7.1, but after an off-chance-grabbing-at-straws upgrade to 2.7.2 ...the answer becomes right! Ugh. Something was amiss with prime()/primepi(), and now apparently fixed.

[CRGreathouse]

Originally Posted by Batalov

It is always interesting when one's answer to a problem is wrong not because the solution is wrong, but because of a bug in Pari.

A certain (rather straightforward) solution for PE 501 gives the wrong answer with Pari 2.7.1, but after an off-chance-grabbing-at-straws upgrade to 2.7.2 ...the answer becomes right! Ugh. Something was amiss with prime()/primepi(), and now apparently fixed.

Can you give more information? I'd like to reproduce this.

[Batalov]

I will PM you. I don't want to "post the solution". (And I have no doubt that you can solve it anyway, and even more elegantly!)

Maybe I'll be able to make the debug case smaller. Anyway - the bug was in 2.7.1, seems to be fixed in 2.7.2. Not the other way around.

EDIT: I removed the code that gives away half of the PE 501 solution, because see below...

[Batalov]

I simplified the bug now to a single call (found in 2.6.0 -- 2.7.1, both linux and windows):

# ./gp -p500000000

? primepi(450450450)
%1 = 23875519  (in 2.7.2, or <=2.5.5)
%1 = 23875520  (in 2.7.1, 2.7.0, 2.6.x with precomputed primes)

# but if you don't precompute primes, the answer is correct, in both cases! 
# the calculation of my script without precomputed primes is too slow, so they have to be precomputed

[ldesnogu]

Originally Posted by Batalov

I simplified the bug now to a single call (found in 2.6.0 -- 2.7.1, both linux and windows):

# ./gp -p500000000

? primepi(450450450)
%1 = 23875519  (in 2.7.2, or <=2.5.5)
%1 = 23875520  (in 2.7.1, 2.7.0, 2.6.x with precomputed primes)

# but if you don't precompute primes, the answer is correct, in both cases! 
# the calculation of my script without precomputed primes is too slow, so they have to be precomputed

Perhaps this?

commit 758f760868e38edb7b4e88f3e1759282623366e6
Author: Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>
Date:   Sat Mar 8 08:41:14 2014 +0100

    29- primepi(N >= 2^32 or 2^64) off by 1
    
    + remove useless HACK

index f622ba2..931abf7 100644
--- a/src/basemath/prime.c
+++ b/src/basemath/prime.c
@@ -1031,22 +1031,24 @@ uprimepi(ulong a)
   {
     byteptr d;
     prime_table_next_p(a, &d, &p, &n);
+    return p == a? n: n-1;
   }
   else
   {
     long i = prime_table_closest_p(a);
     forprime_t S;
     p = prime_table[i].p;
-    if (p > maxp)
+    if (p > a)
     {
       i--;
       p = prime_table[i].p;
     }
+    /* p = largest prime in table <= a */
     n = prime_table[i].n;
     (void)u_forprime_init(&S, p+1, a);
     for (; p; n++) p = u_forprime_next(&S);
+    return n-1;
   }
-  return p == a? n: n-1;
 }
 
 GEN
@@ -1061,7 +1063,6 @@ primepi(GEN x)
   if (signe(N) <= 0) return gen_0;
   avma = av; l = lgefint(N);
   if (l == 3) return utoi(uprimepi(N[2]));
-  new_chunk(l); /*HACK*/
   i = prime_table_len-1;
   p = prime_table[i].p;
   n = prime_table[i].n;
@@ -1069,7 +1070,7 @@ primepi(GEN x)
   nn = setloop(utoipos(n));
   pp = gen_0;
   for (; pp; incloop(nn)) pp = forprime_next(&S);
-  avma = av; return icopy(nn);
+  return gerepileuptoint(av, subiu(nn,1));
 }
 
 /* pi(x) < x/log x * (1 + 1/log x + 2.51/log^2 x)), x>=355991 [ Dusart ]

Edit: That doesn't seem to be it since 2.7.0 was released after that change.

[ldesnogu]

Next try

[ldesnogu@localhost pari]$ git show c2c67fa04290844d8dbd719cdf1d473c4463636c
commit c2c67fa04290844d8dbd719cdf1d473c4463636c
Author: Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>
Date:   Mon May 26 11:56:00 2014 +0200

    25- (gp -p N) or (primelimit=N in gprc_ for N >= 436273290 resulted in an
       incorrect primetable. N.B. Such commands are now useless: needed primes
       are produced dynamically anyway.

diff --git a/CHANGES b/CHANGES
index 551036b..d13ec2f 100644
--- a/CHANGES
+++ b/CHANGES
@@ -7,6 +7,9 @@ Done for version 2.7.2 (released ??/??/2014):
   Fixed
     1- gaffsg(0, t_PADIC): wrong valuation                                [F21]
     2- (t_INTMOD with word-sized modulus)^(huge negative power) [#1584]   [F24]
+    3- (gp -p N) or (primelimit=N in gprc_ for N >= 436273290 resulted in an
+       incorrect primetable. N.B. Such commands are now useless: needed primes
+       are produced dynamically anyway.                                   [F25]
 
 Done for version 2.7.1 (released 16/05/2014):
 [last column crossreferences current development release 2.8.0]
diff --git a/src/basemath/arith2.c b/src/basemath/arith2.c
index bc650ff..5a671c7 100644
--- a/src/basemath/arith2.c
+++ b/src/basemath/arith2.c
@@ -358,7 +358,7 @@ sieve_chunk(byteptr known_primes, ulong s, byteptr data, ulo
   }
 }
 
-/* assume maxnum <= 436273290 < 2^29 */
+/* assume maxnum <= 436273289 < 2^29 */
 static void
 initprimes0(ulong maxnum, long *lenp, ulong *lastp, byteptr p1)
 {
@@ -435,16 +435,18 @@ maxprime(void) { return diffptr ? _maxprime : 0; }
 void
 maxprime_check(ulong c) { if (_maxprime < c) pari_err_MAXPRIME(c); }
 
-/* We ensure 65302 <= maxnum <= 436273290: the LHS ensures modular function
- * have enough fast primes to work, the RHS ensures that p_{n+1} - p_n < 255 */
+/* We ensure 65302 <= maxnum <= 436273289: the LHS ensures modular function
+ * have enough fast primes to work, the RHS ensures that p_{n+1} - p_n < 255
+ * (N.B. RHS would be incorrect since initprimes0 would make it odd, thereby
+ * increasing it by 1) */
 byteptr
 initprimes(ulong maxnum, long *lenp, ulong *lastp)
 {
   byteptr t;
   if (maxnum < 65537)
     maxnum = 65537;
-  else if (maxnum > 436273290)
-    maxnum = 436273290;
+  else if (maxnum > 436273289)
+    maxnum = 436273289;
   t = (byteptr)pari_malloc((size_t) (1.09 * maxnum/log((double)maxnum)) + 146);
   initprimes0(maxnum, lenp, lastp, t);
   return (byteptr)pari_realloc(t, *lenp);