Be Geeks! (Editorial)

The number of pairs $(i, j)$ equals $n \cdot (n + 1) /2$ n \cdot (n + 1) / 2, and for $n = 2 \cdot 1 0^{5}$ n = 2 \cdot 10^5 it is about $2 \cdot 1 0^{10}$ 2 \cdot 10^{10}, so enumerating all segments is impossible. We compute the sum using the "divide and conquer" method, splitting the array at the position of the maximum element.

Splitting at the maximum. Suppose we need to find the sum of $P (i, j)$ P(i, j) over all segments $[i; j]$ lying inside the segment $[b; e]$ . Find the position $m$ of the maximum element on $[b; e]$ . Every segment $[i; j] \subseteq [b; e]$ [i; j] \subseteq [b; e] belongs to exactly one of three types:

it contains position $m$ , that is, $b \leq i \leq m \leq j \leq e$ b \le i \le m \le j \le e;
it lies entirely in $[b; m - 1]$ [b; m - 1];
it lies entirely in $[m + 1; e]$ [m + 1; e].

Segments of the second and third types are processed recursively. For segments of the first type the maximum equals $a_{m}$ a_m, so $M (i, j) = a_{m}$ M(i, j) = a_m, and their total contribution equals

a_{m} \cdot \sum G (i, j), where the sum is taken over all b \leq i \leq m \leq j \leq e

a_m \cdot \sum G(i, j), \quad \text{where the sum is taken over all } b \le i \le m \le j \le e

It remains to compute quickly the sum of GCDs over all segments containing position $m$ .

Decomposing the GCD. A segment $[i; j]$ containing $m$ consists of the left part $[i; m]$ and the right part $[m; j]$ . Therefore

G (i, j) = g cd (G (i, m), G (m, j))

G(i, j) = \gcd\bigl(G(i, m), G(m, j)\bigr)

Thus, the GCD of a segment is determined by a pair of values: the GCD of the left part and the GCD of the right part.

Groups of endpoints. As $j$ increases, the value $G (m, j)$ G(m, j) does not increase. Moreover, whenever it changes, the new value is a proper divisor of the previous one and is therefore at least two times smaller. Hence $G (m, j)$ G(m, j) takes at most $⌊ lo g_{2} C ⌋ + 1$ \lfloor \log_2 C \rfloor + 1 distinct values, where $C$ is the largest element value. For $C \leq 1 0^{9}$ C \le 10^9 this is at most $⌊ lo g_{2} C ⌋ + 1 = 30$ \lfloor \log_2 C \rfloor + 1 = 30 values.

Positions $j$ with the same value of $G (m, j)$ G(m, j) form a contiguous block. We call such a block of positions a group. Each group is described by a pair $(g, cnt)$ (g, \mathit{cnt}), where $g$ is the GCD value and $cnt$ is the number of positions in the group. Similarly, the left endpoints $i$ are split into groups by the value of $G (i, m)$ .

The end of a group is found by binary search. If a group starts at position $j$ and $g = G (m, j)$ g = G(m, j), then the equality $G (m, t) = g$ G(m, t) = g holds for all $t$ from $j$ to the end of the group and fails for all $t$ after it.

Sum over pairs of groups. Let the left group be $(g_{L}, cnt_{L})$ (g_L, \mathit{cnt}_L), and the right group be $(g_{R}, cnt_{R})$ (g_R, \mathit{cnt}_R). Then all $cnt_{L} \cdot cnt_{R}$ \mathit{cnt}_L \cdot \mathit{cnt}_R segments whose left endpoint lies in the left group and whose right endpoint lies in the right group have the same GCD $g cd (g_{L}, g_{R})$ \gcd(g_L, g_R). Therefore

\sum G (i, j) = \sum g cd (g_{L}, g_{R}) \cdot cnt_{L} \cdot cnt_{R},

\sum G(i, j) = \sum \gcd(g_L, g_R) \cdot \mathit{cnt}_L \cdot \mathit{cnt}_R,

where the second sum is taken over all pairs of a left and a right group. The number of such pairs does not exceed $30 \cdot 30$ 30 \cdot 30.

Range queries. To find the position of the maximum and to compute the GCD on a segment, we use two sparse tables. Both operations are idempotent: $max (x, x) = x$ \max(x, x) = x and $g cd (x, x) = x$ \gcd(x, x) = x. Therefore a range query is answered using two overlapping blocks of length $2^{k}$ 2^k.

The value of $F (A)$ is very large, so all computations are performed modulo $1 0^{9} + 7$ 10^9 + 7.

Complexity. Each recursive call removes one position $m$ from consideration, so there are at most $n$ calls. In each call there are at most $O (lo g C)$ O(\log C) groups on each side. The end of a group is found by binary search with $O (lo g n)$ O(\log n) GCD queries, and each query takes $O (lo g C)$ O(\log C). Iterating over pairs of groups requires $O (lo g^{2} C)$ O(\log^2 C) GCD computations. Building the sparse tables takes $O (n lo g n \cdot lo g C)$ O(n \log n \cdot \log C). The overall complexity is $O (n lo g^{2} C \cdot (lo g n + lo g C))$ O(n \log^2 C \cdot (\log n + \log C)). In practice there are far fewer than $30$ groups, and the algorithm runs fast.

Example

Consider the second example: $A = (2, 4, 6, 12, 3)$ A = (2, 4, 6, 12, 3). Positions are numbered from one.

The first call processes the whole array $[1; 5]$ . The maximum element is $a_{4} = 12$ a_4 = 12, so $m = 4$ m = 4. The figure shows the groups of left and right endpoints.

Left endpoints: $G (4, 4) = 12$ G(4, 4) = 12, $G (3, 4) = g cd (6, 12) = 6$ G(3, 4) = \gcd(6, 12) = 6, $G (2, 4) = G (1, 4) = 2$ G(2, 4) = G(1, 4) = 2. We get the groups $(12, 1)$ , $(6, 1)$ and $(2, 2)$ .

Right endpoints: $G (4, 4) = 12$ G(4, 4) = 12, $G (4, 5) = g cd (12, 3) = 3$ G(4, 5) = \gcd(12, 3) = 3. We get the groups $(12, 1)$ and $(3, 1)$ .

Let us go through all pairs of groups:

left group	right group	GCD	segments	contribution
(12, 1)	(12, 1)	12	1	12
(12, 1)	(3, 1)	3	1	3
(6, 1)	(12, 1)	6	1	6
(6, 1)	(3, 1)	3	1	3
(2, 2)	(12, 1)	2	2	4
(2, 2)	(3, 1)	1	2	2

The sum of GCDs over the segments containing position $4$ is $12 + 3 + 6 + 3 + 4 + 2 = 30$ 12 + 3 + 6 + 3 + 4 + 2 = 30, and the contribution of these segments is $12 \cdot 30 = 360$ 12 \cdot 30 = 360. For example, the last row of the table corresponds to the segments $[1; 5]$ and $[2; 5]$ : for each of them the GCD is 1 and the maximum is 12.

Then the segments $[1; 3]$ and $[5; 5]$ are processed recursively. All calls are listed in the table:

segment	$m$	$a_{m}$ \boldsymbol{a_m}	left groups	right groups	GCD sum	contribution
$[1; 5]$ \boldsymbol{[1; 5]}	4	12	(12, 1), (6, 1), (2, 2)	(12, 1), (3, 1)	30	360
$[1; 3]$ \boldsymbol{[1; 3]}	3	6	(6, 1), (2, 2)	(6, 1)	10	60
$[1; 2]$ \boldsymbol{[1; 2]}	2	4	(4, 1), (2, 1)	(4, 1)	6	24
$[1; 1]$ \boldsymbol{[1; 1]}	1	2	(2, 1)	(2, 1)	2	4
$[5; 5]$ \boldsymbol{[5; 5]}	5	3	(3, 1)	(3, 1)	3	9

The answer is $360 + 60 + 24 + 4 + 9 = 457$ 360 + 60 + 24 + 4 + 9 = 457.

Algorithm Implementation

Declare the constants:

$MOD = 1 0^{9} + 7$ \mathit{MOD} = 10^9 + 7 — the modulus,
$LOG = 18$ \mathit{LOG} = 18 — the largest exponent in the sparse tables ( $2^{18} > 2 \cdot 1 0^{5}$ 2^{18} > 2 \cdot 10^5).

Declare the type ll.

#define MOD 1000000007LL
#define LOG 18 // 2^18 = 262144 > 2*10^5
typedef long long ll;

The sequence is stored in the array $a$ ; elements are indexed from zero. The value $lg2 [len]$ equals $⌊ lo g_{2} len ⌋$ \lfloor \log_2 \mathit{len} \rfloor. The sparse table $spMax [i] [k]$ \mathit{spMax}[i][k] stores the index of the maximum element, and $spGcd [i] [k]$ \mathit{spGcd}[i][k] stores the GCD of the elements on the segment $[i; i + 2^{k} - 1]$ [i; i + 2^k - 1].

vector<int> a;
vector<int> lg2;          // lg2[len] = floor(log2(len))
// spMax[i][k] — index of the maximum on [i, i + 2^k)
vector<vector<int>> spMax;
// spGcd[i][k] — GCD on [i, i + 2^k)
vector<vector<int>> spGcd;

The function gcd computes the GCD of two numbers using the Euclidean algorithm.

int gcd(int x, int y)
{
  while (y) { int t = x % y; x = y; y = t; }
  return x;
}

The function better, given two indices $i$ and $j$ , returns the one whose array element is larger. If the elements are equal, it returns $i$ .

int better(int i, int j)
{
  return a[i] >= a[j] ? i : j;
}

The function build fills the array $lg2$ and builds both sparse tables in $O (n lo g n \cdot lo g C)$ O(n \log n \cdot \log C), $C \leq 1 0^{9}$ C \le 10^9.

void build()
{
  lg2.assign(n + 1, 0);
  for (int len = 2; len <= n; len++) lg2[len] = lg2[len / 2] + 1;
  spMax.assign(n, vector<int>(LOG + 1));
  spGcd.assign(n, vector<int>(LOG + 1));
  for (int i = 0; i < n; i++)
  {
    spMax[i][0] = i;
    spGcd[i][0] = a[i];
  }
  for (int k = 1; (1 << k) <= n; k++)
    for (int i = 0; i + (1 << k) <= n; i++)
    {
      int h = 1 << (k - 1);    // half of the block length
      spMax[i][k] = better(spMax[i][k - 1], spMax[i + h][k - 1]);
      spGcd[i][k] = gcd(spGcd[i][k - 1], spGcd[i + h][k - 1]);
    }
}

C++

The function qMax returns the index of the maximum element on the segment $[l; r]$ [l; r].

int qMax(int l, int r)
{
  int k = lg2[r - l + 1];
  return better(spMax[l][k], spMax[r - (1 << k) + 1][k]);
}

The function qGcd returns the GCD of the elements on the segment $[l; r]$ [l; r]. Blocks of length $2^{k}$ 2^k may overlap, but this does not affect the result, since $g cd (x, x) = x$ \gcd(x, x) = x.

int qGcd(int l, int r)
{
  int k = lg2[r - l + 1];
  return gcd(spGcd[l][k], spGcd[r - (1 << k) + 1][k]);
}

The function lastRight finds the rightmost position $t \in [from; e]$ t \in [\mathit{from}; e] such that $G (m, t) = g$ G(m, t) = g. Here $g = G (m, from)$ g = G(m, \mathit{from}). Since $G (m, t)$ does not increase as $t$ grows, the equality holds on some initial part of the segment $[from; e]$ , and binary search finds its last position. The midpoint is computed as $(lo + hi + 1) /2$ (\mathit{lo} + \mathit{hi} + 1) / 2 so that the search does not get stuck on the assignment $lo = mid$ \mathit{lo} = \mathit{mid}.

int lastRight(int m, int from, int e, int g)
{
  int lo = from, hi = e;
  while (lo < hi)
  {
    int mid = (lo + hi + 1) / 2;
    if (qGcd(m, mid) == g) lo = mid; else hi = mid - 1;
  }
  return lo;
}

C++

The function firstLeft symmetrically finds the leftmost position $t \in [b; from]$ t \in [b; \mathit{from}] such that $G (t, m) = g$ G(t, m) = g.

int firstLeft(int m, int from, int b, int g)
{
  int lo = b, hi = from;
  while (lo < hi)
  {
    int mid = (lo + hi) / 2;
    if (qGcd(mid, m) == g) hi = mid; else lo = mid + 1;
  }
  return lo;
}

C++

The structure Group describes a group of segment endpoints:

$g$ — the GCD value,
$cnt$ — the number of positions in the group.

struct Group { int g; ll cnt; };

The function solve returns, modulo $MOD$ , the sum of $P (i, j)$ P(i, j) over all segments $[i; j]$ lying inside $[b; e]$ . For an empty segment the sum is $0$ .

ll solve(int b, int e)
{
  if (b > e) return 0;

Find the position $m$ of the maximum element on $[b; e]$ .

  int m = qMax(b, e);

Split the right endpoints $j = m, m + 1, \dots, e$ j = m, m + 1, \ldots, e into groups. For the current $j$ compute $g = G (m, j)$ g = G(m, j) and use the function lastRight to find the last position $t$ of the group. The group contains $t - j + 1$ t - j + 1 positions; the next group starts at position $t + 1$ t + 1.

  vector<Group> R;
  for (int j = m; j <= e; )
  {
    int g = qGcd(m, j);
    int t = lastRight(m, j, e, g);
    R.push_back({ g, t - j + 1 });
    j = t + 1;
  }

C++

Similarly, split the left endpoints $i = m, m - 1, \dots, b$ i = m, m - 1, \ldots, b into groups. The group occupies positions from $t$ to $i$ and contains $i - t + 1$ i - t + 1 positions; the next group ends at position $t - 1$ t - 1.

  vector<Group> L;
  for (int i = m; i >= b; )
  {
    int g = qGcd(i, m);
    int t = firstLeft(m, i, b, g);
    L.push_back({ g, i - t + 1 });
    i = t - 1;
  }

C++

Compute the sum $sumG$ of GCD values over all segments containing position $m$ . For each pair of groups add $g cd (l . g, r . g) \cdot l . cnt \cdot r . cnt$ \gcd(l.g, r.g) \cdot l.\mathit{cnt} \cdot r.\mathit{cnt}. The product $l . cnt \cdot r . cnt$ l.\mathit{cnt} \cdot r.\mathit{cnt} does not exceed $4 \cdot 1 0^{10}$ 4 \cdot 10^{10} and fits into the type long long. We reduce it modulo $MOD$ ; after that, the product of the GCD (less than $2^{30}$ 2^{30}) and a number less than $MOD$ also fits into long long.

  ll sumG = 0;
  for (Group &l : L)
    for (Group &r : R)
    {
      ll cnt = (l.cnt * r.cnt) % MOD;
      sumG = (sumG + gcd(l.g, r.g) % MOD * cnt) % MOD;
    }

C++

The maximum on all segments containing $m$ equals $a_{m}$ a_m, so their total contribution equals $a_{m} \cdot sumG$ a_m \cdot \mathit{sumG}.

  ll res = (ll)a[m] % MOD * sumG % MOD;

Add the recursively computed sums for the segments lying entirely to the left and to the right of position $m$ .

  return (res + solve(b, m - 1) + solve(m + 1, e)) % MOD;
}

The main part of the program. Read the input data.

scanf("%d", &n);
a.resize(n);
for (int &x : a) scanf("%d", &x);

Build the sparse tables.

build();

Print the answer for the whole array — the segment $[0; n - 1]$ [0; n - 1].

printf("%lld\n", solve(0, n - 1));