題目

Problem 585: Nested square roots

Consider the term that is representing a nested square root. x, y and z are positive integers and y and z are not allowed to be perfect squares, so the number below the outer square root is irrational. Still it can be shown that for some combinations of x, y and z the given term can be simplified into a sum and/or difference of simple square roots of integers, actually denesting the square roots in the initial expression.

Here are some examples of this denesting:

As you can see the integers used in the denested expression may also be perfect squares resulting in further simplification.

Let F(n) be the number of different terms , that can be denested into the sum and/or difference of a finite number of square roots, given the additional condition that 0 < x ≤ n. That is,

with k, x, y, z and all a_i being positive integers, all s_i = ±1 and x ≤ n.
Furthermore y and z are not allowed to be perfect squares.

Nested roots with the same value are not considered different, for example , and , that can all three be denested into , would only be counted once.

You are given that F(10) = 17, F(15) = 46, F(20) = 86, F(30) = 213 and F(100) = 2918 and F(5000) = 11134074.
Find F(5000000).

分析

我們有以下情形（以下各變數都是正整數）：

1. 
   當 i = 1, b = 2 時，就是題目中的第 1 個例子。
2. 
   當 a = 3, b = 5 時，就是題目中的第 2 個例子。
   Update: 當 a = 8, b = 18 時，這是有問題的：
   
3. 
4. 
   Update: i = 2, c = 6, a = 2.5 也是可以的。
5. 
   當 a = 3, b = 2, c = 3 時，就是題目中的第 3 個例子。
   當 a = 3, b = 2, c = 5 時，就是題目中的第 4 個例子。
   Update: a = 1.5, b = 6, c = 8 也是可以的。

Update: 使用 SymPy，得到以下幾個不屬於前面 5 種情形的例子：

在情形05中，允許 a = b 或者 a = c。
在情形05中，如果令 b = i²，就變為情形04，此時，不允許 a = c。
在情形04中，如果令 i = 1，就變為情形03，此時，要求 c > a。
在情形02中，如果令 a = i²，就變為情形01，此時，不要求 b > i²。

不成功的解答

根據以上分析，我們有以下 C 程式：

 1: #include <stdio.h>
 2: #include <math.h>
 3: 
 4: int isSquare(int n) { double x = sqrt(n); return x == (int)x; }
 5: int max(int a, int b) { return (a > b) ? a : b; }
 6: 
 7: int count(int b, int c0, int c9)
 8: { // count: c in (c0, c9]; c,c/b != k^2
 9:   return c9 - (int)sqrt(c9) - (int)sqrt(c9/b)
10:        - c0 + (int)sqrt(c0) + (int)sqrt(c0/b);
11: }
12: 
13: int count4(int a, int b, int c0, int c9)
14: { // count: c in (c0, c9]; c,c/b,ac/b != k^2
15:   int z = 0;
16:   for (int d, c = c0 + 1; c <= c9; c++) {
17:     if (isSquare(c)) continue;
18:     if (c % b == 0 && isSquare(c / b)) continue;
19:     if ((d = a * c) % b == 0 && isSquare(d / b)) continue;
20:     z++;
21:   }
22:   return z;
23: }
24: 
25: long f2a(int n)
26: { // i + √b: 0 < i, 1 < b != k^2
27:   long z = 0;    // i^2 + b <= n
28:   for (int b, i = 1; 1 < (b = n - i * i); i++)
29:     z += b - (int)sqrt(b);
30:   return z;
31: }
32: 
33: long f2b(int n)
34: { // √a + √b: 1 < a < b; a,b,b/a != k^2
35:   long z = 0;             // a + b <= n
36:   for (int b, a = 2; a < (b = n - a); a++)
37:     if (!isSquare(a)) z += count(a, a, b);
38:   return z;
39: }
40: 
41: long f4a(int n)
42: { // -1 + √c + √a + √(ac): c > a > 1; a,c,c/a != k^2
43:   long z = 0;                 // (a + 1)(c + 1) <= n
44:   for (int c, a = 2; a < (c = n / (a + 1) - 1); a++)
45:     if (!isSquare(a)) z += count(a, a, c);
46:   return z;
47: }
48: 
49: long f4b(int n)
50: { // -i + √c + i√a + √(ac):   c > i^2 > 1; a,c != k^2
51:   long z = 0; // a != c, a > 1, (a + 1)(c + i^2) <= n
52:   for (int i2, i = 2; (i2 = i * i) * 2 * 3 <= n; i++) {
53:     for (int b, c, a = 2; (b = max(a,i2)) < (c = n/(a+1)-i2); a++)
54:       if (!isSquare(a)) z += count(a, b, c); // c > a; c/a != k^2
55:     for (int a, c = i2 + 1; c < (a = n / (i2 + c) - 1); c++)
56:       if (!isSquare(c)) z += count(c, c, a); // a > c; a/c != k^2
57:   }
58:   return z;
59: }
60: 
61: long f4c(int n)
62: { // -√b + √c + √(ab) + √(ac):   a,b,c,c/b,ac/b != k^2
63:   long z = 0; // a > 1, c > b > 1, (a + 1)(b + c) <= n
64:   for (int a = 2; 5 * (a + 1) <= n; a++)
65:     if (!isSquare(a)) // a == (b or c) is allowed
66:       for (int c, b = 2; b < (c = n / (a + 1) - b); b++)
67:         if (!isSquare(b)) z += count4(a, b, b, c);
68:   return z;
69: }
70: 
71: long f(int n)
72: {
73:   return f2a(n) + f2b(n) + f4a(n) + f4b(n) + f4c(n);
74: }
75: 
76: int main()
77: {
78:   static int a[] = { 10, 15, 20, 30, 100, 5000 };
79:   for (int i = 0; i < sizeof(a) / sizeof(a[0]); i++)
80:     printf("%ld ", f(a[i])); // 2813  8237158
81:   printf(       "\n17 46 86 213 2918 11134074\n");
82: }

這個程式執行結果如下：

17 46 86 213 2813 8237158
17 46 86 213 2918 11134074

前 4 個值符合，後 2 個值不符。肯定是漏掉某些情形。

可以工作的程式

Update: 經過這幾天的思考，得到以下 C 程式：

 1: #include <stdio.h>
 2: #include <math.h>
 3: 
 4: int isSquare(long n) { long x = (long)sqrt(n); return x * x == n; }
 5: int gcd(int a, int b) { return b == 0 ? a : gcd(b, a % b); }
 6: 
 7: int count(int a, int b0, int b9)
 8: { // count: b in (b0, b9]; b,ab != k^2
 9:   for (int k, j = a % 2 + 1, i = j + 1; (k = i * i) <= a; )
10:     if (a % k == 0) a /= k; else i += j;
11:   return b9 - (int)sqrt(b9) - (int)sqrt(b9/a)
12:        - b0 + (int)sqrt(b0) + (int)sqrt(b0/a);
13: }
14: 
15: long f2a(int n)
16: { // i + √b: 0 < i, 1 < b != k^2
17:   long z = 0;    // i^2 + b <= n
18:   for (int b, i = 1; 1 < (b = n - i * i); i++)
19:     z += b - (int)sqrt(b);
20:   return z;
21: }
22: 
23: long f2b(int n)
24: { // √a + √b: 1 < a < b; a,b,ab != k^2
25:   long z = 0;            // a + b <= n
26:   for (int b, a = 2; a < (b = n - a); a++)
27:     if (!isSquare(a)) z += count(a, a, b);
28:   return z;
29: }
30: 
31: long f4(int n)
32: { // -√a + √b + √c + √d: ad = bc;     ab,ac,bc != k^2
33:   long z = 0;    // 0 < a < b < c, a + b + c + d <= n
34:   for (int a = 1; (a << 2) + 6 <= n; a++)
35:     for (int b = a + 1; b < (long)a*(n-a-b)/(a+b); b++) {
36:       if (isSquare((long)a * b)) continue;
37:       int k = a / gcd(a, b), c = (b + 1) / k * k;
38:       if (c < b + 1) c += k;
39:       for (long e; a + b + c + c + 1 <= n; c += k)
40:         if (a + b + c + (e = (long)c * b) / a <= n &&
41:           !isSquare(e) && !isSquare((long)c * a)) z++;
42:     }
43:   return z;
44: }
45: 
46: long f(int n) { return f2a(n) + f2b(n) + f4(n); }
47: 
48: int main()
49: {
50:   static int a[] = { 10, 15, 20, 30, 100, 5000 };
51:   for (int i = 0; i < sizeof(a) / sizeof(a[0]); i++)
52:     printf("%ld ", f(a[i]));
53:   printf("\n17 46 86 213 2918 11134074\n");
54: }

這個程式的執行時間是 0.3 秒，執行結果是：

17 46 86 213 2918 11134074
17 46 86 213 2918 11134074

這個結果與題目中給出的值一致，應該是正確的。但是用來計算 F(5000000) 就太慢了。

在一臺 16 個邏輯 CPU 的 openSUSE Tumbleweed 64-bit 的伺服器上，這個程式稍做修改，使用 OpenMP 進行平行計算，執行了約 9 個小時，終於得到正確的結果。