euler/src/p44.cpp at master · fancidev/euler · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
/**
 * Pentagonal numbers are generated by the formula, P(n)=n(3n-1)/2.
 * The first ten pentagonal numbers are:
 *
 *   1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ...
 *
 * It can be seen that P(4) + P(7) = 22 + 70 = 92 = P(8). However,
 * their difference, 70 - 22 = 48, is not pentagonal.
 *
 * Find the pair of pentagonal numbers, P(j) and P(k), for which
 * their sum and difference is pentagonal and D = |P(k) - P(j)| is
 * minimised. What is the value of D?
 *
 * SOLUTION
 *
 * This program finds the following D(r) for each r-gonal numbers
 * for r = 3 to 16, except for r = 4. Finding D(11) takes the most
 * time (around 9 seconds).
 *
 * D(3)  = 6
 * D(5)  = 5482660
 * D(6)  = 1633528
 * D(7)  = 107017308
 * D(8)  = 1281
 * D(9)  = 2197072950
 * D(10) = 110360277
 * D(11) = 8864150212750
 * D(12) = 177048201
 * D(13) = 1023061390
 * D(14) = 695515899
 * D(15) = 30653498337
 * D(16) = 46941160
 */

#define EULER_LOOPLESS_DIVISOR_GENERATION 1
#include <cstdint>
#include <iostream>
#include "euler/prime_factor.hpp"
#include "euler/divisor.hpp"
#include "euler/figurate.hpp"
#include "euler.h"

BEGIN_PROBLEM(44, solve_problem_44)
  PROBLEM_TITLE("Pentagonal numbers whose sum and difference are also pentagonal")
  PROBLEM_ANSWER("5482660")
  PROBLEM_DIFFICULTY(2)
  PROBLEM_FUN_LEVEL(3)
  PROBLEM_TIME_COMPLEXITY("n^2")
  PROBLEM_SPACE_COMPLEXITY("log(n)")
  PROBLEM_KEYWORDS("polygonal number")
END_PROBLEM()

// Find the minimum difference between two r-gonal numbers
// whose sum and difference are both r-gonal numbers.
static int64_t find_min_diff(int r)
{
  int count1 = 0, count2 = 0;

  typedef int64_t calc_t;
  for (int i = 1; ; i++)
  {
    calc_t D = static_cast<calc_t>(i) * static_cast<calc_t>((r-2)*i-(r-4));

    // Enumerate the divisors of D.
    auto f = euler::merge(
        euler::factorize(i), euler::factorize((r-2)*i-(r-4)));
    for (calc_t a: euler::divisors<calc_t>(f))
    {
      ++count1;
      calc_t b = D/a;
      if (a >= b)
      {
        continue;
      }
      if ((b+(r-4))%(r-2) != 0)
      {
        continue;
      }
      b=(b+(r-4))/(r-2);
      if ((a+b)%2 != 0)
      {
        continue;
      }
      calc_t j = (a+b)/2;
      calc_t k = (b-a)/2;
      if (k <= 0)
      {
        continue;
      }

      // Check whether P_j + P_k is a polygonal number of rank r.
      ++count2;
      calc_t Pl = euler::polygonal(r, j) + euler::polygonal(r, k);
      if (euler::is_polygonal(r, Pl))
      {
        if (verbose())
        {
          std::cout << "D(" << r << ") = " << D/2 << ", i,j,k = "
            << i << "," << j << "," << k << std::endl;
          std::cout << "    # divisors checked: " << count1 << std::endl;
          std::cout << "    # polygonal tested: " << count2 << std::endl;
        }
        return D/2;
      }
    }
  }
}

static void solve_problem_44()
{
#if 0
  for (int r = 3; r <= 16; r++)
  {
    if (r != 4)
    {
      find_min_diff(r);
    }
  }
#else
  std::cout << find_min_diff(5) << std::endl;
#endif
}