PROBLEM LINK:
Practice
Contest: Division 1
Contest: Division 2
Setter: Aleksa Plasvic
Tester: Hasan Jaddouh
Editorialist: Taranpreet Singh
DIFFICULTY:
Hard
PREREQUISITES:
Number theory, partial sums, Observations and Branch and bound.
PROBLEM:
Given a sequence A of length N. GCD value of the sequence is the number of non-empty subarrays which have GCD of all values greater than one. We want to maximize the GCD value of the sequence by changing at most one element to any value in the range [1, 5*10^5].
QUICK EXPLANATION
- Any number up to 5*10^5 can have at most 6 distinct prime factors since the product of first seven primes exceed 5*10^5.
- For any position, we try to replace it to any value which has prime factors same as the two adjacent elements. So, we try to select the subset of primes out of 11 distinct primes such that their product does not exceed 5*10^5 and the number of subarrays is maximum.
- For calculating the number of subarrays, we can calculate it in three parts, Number of subarrays ending before position p, Number of subarrays starting after position p and the number of subarrays covering position p. Count of subarrays of the first two types can be calculated pre-calculated.
- For calculating subarrays of the third type, considering each subset of primes out of 11 primes, such that their product is up to 5*10^5, then for each prime, we find the range of elements containing position p all of which have that prime as their factor. We get at most 11 such ranges. We need to count the number of distinct subarrays out of these ranges which contain the position p and lie within at least one range. Branch and bound can be used to speed up the selection of valid subsets.
EXPLANATION
For each position, we try to replace that value with the number of subarrays and try to find the maximum GCD value of an array.
For each position p, we can divide the number of subarrays into three parts, subarrays ending before position p, subarrays starting after position p and subarrays containing position p. It can be seen that the first two categories aren’t affected by value at position p and can be precalculated.
See, Let us find the number of good subarrays ending at position p. Suppose for any prime factor of A_p, if x immediately before p elements including A_p have p as their prime factor, we can choose any of the x as starting point of the subarray. Now, the number of subarrays ending before or at position p is just the prefix sum array of these values. We can do the same for subarrays starting after any position. For finding x for each position, we can make prefix and suffix maps for each position, where prefix maps store the prime factors of current position as key, and number of elements immediately to the left including that element, which has that prime as their factor.
Now, calculating subarrays of the third type require some effort.
It is easy to notice that if we can choose ith value as the product (Or LCM) of two adjacent values, the number of subarrays of the third type shall be maximized. But we cannot do so since we need the value to be in the range [1, 5*10^5] and LCM may go up to 25*10^{10} in the worst case.
So, we have to choose a subset of primes such that their product does not exceed 5*10^5. Also, it doesn’t make sense to include any prime multiple times, as we do not care about actual GCD of subarray as long as it is greater than one.
Now, for each prime, we can find the length of the segment to the immediate left and right which have that prime as their factor. Now, all subarrays which lie within the range for at least one of the prime chosen is a good subarray of type three.
Consider example
Given array 4 6 1 15 and we are checking all subsets of primes for position 3. 6 has factors 2 and 3 while 15 has factors 3 and 5. So, we try all subsets of setting (2,3,5) and see the best subset.
We see, that prime two appear in the range [1, 3], prime 3 appears in the range [2, 4] and prime 5 appears in the range [3,4].
So suppose we are considering subset (2,3,5). Good subarrays which contain position 3 and are included in at least one range are (1,3), (2,3), (2,4), (3,3), (3,4). we need a way to count these subarrays fast for each subset of primes and if the product of selected primes is up to max value, we find the maximum number of subarrays and print answer.
So, now we have at max 11 intervals and a position p, we want to count the number of distinct subarrays which are covered at least one of the selected interval.
To avoid overcounting any subarray, we sort the intervals by the left end and keep a value prev (initially assigned p-1) which says that all intervals ending before or at position prev are already considered. Now, considering an interval [L_i, R_i], and if R_i > prev, there are R_i - prev positions which can be chosen as right end of subarray. Now, for every such right end, we can choose left end any value in the range [L_i, p] since all of these would not have been already counted and all includes p. So, Number of good subarrays found till right end R_i is (i-L_i+1)*(R_i-prev). Now we have considered subarrays ending till R_i, so prev gets updated to R_i.
Repeating this process allows us to calculate the number of subarrays which are contained in at least one range.
Code for Number of distinct subarrays
prev = p-1;
ans = 0;
for interval (L, R): sorted(allIntervals)
if R > prev
ans += (R-prev)*(p-L+1)
prev = R
The same trick was used in problem YVSTR
Now, Let us analyze complexity. For every N positions, we have 6 primes, and then for each position, we consider 2^{11} subsets and need to iterate over 11 intervals for each subset, plus the log factors from prefix and suffix maps, Leading to complexity O(N*C*2^C*log(C)) with a large constant factors due to maps as well as some other where C = 11, which is too much.
Here Setter’s approach is to use the branch and bound as there are many subsets which have product greater than the max value (All subsets with 7 or more primes are automatically invalid, at the very least). So, we recurse, building only valid subsets and calculating the number of subarrays for those cases only. It was found that there were only around 500 such combinations in the worst case.
Editorialist, who’s determined to try all subsets, removed the log factors and a lot of complicated optimizations (took well over three hours, so not recommended).
Time Complexity
Time complexity is O(N*6*f(C)*C*log(C)) where C = 11, the maximum number of factors any value up to 25*10^{10} can have and f(C) denote the number of combinations in the setter’s solution.
SOLUTIONS:
Setter's Solution
#include<bits/stdc++.h>
using namespace std;
#define pb push_back
const int maxN = 5e5+5;
const int maxA = 500000;
struct interval
{
int l;
int r;
};
struct par
{
interval i;
int p;
};
int n;
vector<int> primes[maxN];
long long beginAns,ans;
int l[6][maxN], r[6][maxN];
int ln[maxN], cr[maxN], bg[maxN], a[maxN], ri[maxN];
interval in[20];
par pi[maxN][12];
int included;
bool cmp(par f, par s)
{
if (f.i.l<s.i.l) return true;
if (f.i.l==s.i.l){
return f.i.l>s.i.l;
}
return false;
}
void factorisation()
{
for (int i=2; i<maxN; i++)
if (!primes[i].size())
{
for (int j=i; j<maxN; j+=i)
primes[j].pb(i);
}
}
long long calcIntervals(interval * in, int mid, int c)
{
long long ans = 0;
c--;
mid--;
int mx = 0;
for (int i=0;i<c;i++)
{
mx=max(mx, in[i].r-mid);
ans=ans+(in[i+1].l-in[i].l)*mx;
}
mx = max(mx, in[c].r-mid);
ans+=(mid-in[c].l+2)*mx;
return ans;
}
void preCalculation()
{
for (int i=1; i<=n; i++)
{
if (a[i]==1)
continue;
int mi = i;
for (int j:primes[a[i]])
{
if (i>1 && a[i-1]%j==0)
mi=min(mi,bg[j]);
else
bg[j]=i;
l[ln[i]][i]=bg[j];
ln[i]++;
}
beginAns+=i-mi+1;
}
for (int i=n; i>0; i--)
{
ln[i]=0;
if (a[i]==1)
continue;
for (int j:primes[a[i]])
{
if (i==n || a[i+1]%j)
bg[j]=i;
r[ln[i]][i]=bg[j];
ln[i]++;
}
}
for (int i=1;i<=n;i++)
{
for (int j=0;j<ln[i-1];j++)
if (a[i+1]%primes[a[i-1]][j]==0)
{
interval inter = {l[j][i-1],0};
for (int k=0;k<ln[i+1];k++)
if (primes[a[i+1]][k]==primes[a[i-1]][j])
inter.r = r[k][i+1];
pi[i][ri[i]]={inter, primes[a[i-1]][j]};
ri[i]++;
} else
{
interval inter = {l[j][i-1], i};
pi[i][ri[i]]={inter, primes[a[i-1]][j]};
ri[i]++;
}
for (int j=0;j<ln[i+1];j++)
if (a[i-1]%primes[a[i+1]][j])
{
interval inter = {i, r[j][i+1]};
pi[i][ri[i]]={inter, primes[a[i+1]][j]};
ri[i]++;
}
sort(pi[i],pi[i]+ri[i], cmp);
int cur = 0;
for (int j=0;j<ri[i];j++)
if (a[i]%pi[i][j].p==0)
in[cur++]=pi[i][j].i;
if (!cur && a[i]>1) cr[i]=1; else
if (a[i]>1)
cr[i]=calcIntervals(in, i, cur);
}
return;
}
void rek(long long x, int cur, int m)
{
if (x>maxA) return;
if (cur==ri[m] || 1ll*pi[m][cur].p*x>maxA){
ans = max(ans, beginAns+calcIntervals(in, m, included)-cr[m]);
return;
}
for (int i=cur;i<ri[m];i++){
in[included++] = pi[m][i].i;
rek(1ll*x*pi[m][i].p , i+1, m);
included--;
}
}
int main()
{
cin>>n;
assert(n>0 && n<=50000);
for (int i=1; i<=n; i++)
{
scanf("%d",&a[i]);
assert(a[i]<=500000 && a[i]>0);
}
a[0] = 1;
a[n+1] = 1;
factorisation();
preCalculation();
for (int i=1;i<=n;i++){
included = 0;
rek(1, 0, i);
}
cout<<ans<<endl;
return 0;
}
Tester's Solution
#include<bits/stdc++.h>
using namespace std;
#define pb push_back
const int maxN = 5e5+5;
const int maxA = 500000;
struct interval
{
int l;
int r;
};
struct par
{
interval i;
int p;
};
int n;
vector<int> primes[maxN];
long long beginAns,ans;
int l[6][maxN], r[6][maxN];
int ln[maxN], cr[maxN], bg[maxN], a[maxN], ri[maxN];
interval in[20];
par pi[maxN][12];
int included;
bool cmp(par f, par s)
{
if (f.i.l<s.i.l) return true;
if (f.i.l==s.i.l){
return f.i.l>s.i.l;
}
return false;
}
void factorisation()
{
for (int i=2; i<maxN; i++)
if (!primes[i].size())
{
for (int j=i; j<maxN; j+=i)
primes[j].pb(i);
}
}
long long calcIntervals(interval * in, int mid, int c)
{
long long ans = 0;
c--;
mid--;
int mx = 0;
for (int i=0;i<c;i++)
{
mx=max(mx, in[i].r-mid);
ans=ans+(in[i+1].l-in[i].l)*mx;
}
mx = max(mx, in[c].r-mid);
ans+=(mid-in[c].l+2)*mx;
return ans;
}
void preCalculation()
{
for (int i=1; i<=n; i++)
{
if (a[i]==1)
continue;
int mi = i;
for (int j:primes[a[i]])
{
if (i>1 && a[i-1]%j==0)
mi=min(mi,bg[j]);
else
bg[j]=i;
l[ln[i]][i]=bg[j];
ln[i]++;
}
beginAns+=i-mi+1;
}
for (int i=n; i>0; i--)
{
ln[i]=0;
if (a[i]==1)
continue;
for (int j:primes[a[i]])
{
if (i==n || a[i+1]%j)
bg[j]=i;
r[ln[i]][i]=bg[j];
ln[i]++;
}
}
for (int i=1;i<=n;i++)
{
for (int j=0;j<ln[i-1];j++)
if (a[i+1]%primes[a[i-1]][j]==0)
{
interval inter = {l[j][i-1],0};
for (int k=0;k<ln[i+1];k++)
if (primes[a[i+1]][k]==primes[a[i-1]][j])
inter.r = r[k][i+1];
pi[i][ri[i]]={inter, primes[a[i-1]][j]};
ri[i]++;
} else
{
interval inter = {l[j][i-1], i};
pi[i][ri[i]]={inter, primes[a[i-1]][j]};
ri[i]++;
}
for (int j=0;j<ln[i+1];j++)
if (a[i-1]%primes[a[i+1]][j])
{
interval inter = {i, r[j][i+1]};
pi[i][ri[i]]={inter, primes[a[i+1]][j]};
ri[i]++;
}
sort(pi[i],pi[i]+ri[i], cmp);
int cur = 0;
for (int j=0;j<ri[i];j++)
if (a[i]%pi[i][j].p==0)
in[cur++]=pi[i][j].i;
if (!cur && a[i]>1) cr[i]=1; else
if (a[i]>1)
cr[i]=calcIntervals(in, i, cur);
}
return;
}
void rek(long long x, int cur, int m)
{
if (x>maxA) return;
if (cur==ri[m] || 1ll*pi[m][cur].p*x>maxA){
ans = max(ans, beginAns+calcIntervals(in, m, included)-cr[m]);
return;
}
for (int i=cur;i<ri[m];i++){
in[included++] = pi[m][i].i;
rek(1ll*x*pi[m][i].p , i+1, m);
included--;
}
}
int main()
{
cin>>n;
assert(n>0 && n<=50000);
for (int i=1; i<=n; i++)
{
scanf("%d",&a[i]);
assert(a[i]<=500000 && a[i]>0);
}
a[0] = 1;
a[n+1] = 1;
factorisation();
preCalculation();
for (int i=1;i<=n;i++){
included = 0;
rek(1, 0, i);
}
cout<<ans<<endl;
return 0;
}
Editorialist's Solution
import java.util.*;
import java.io.*;
import java.text.*;
//Solution Credits: Taranpreet Singh
class GCDS{
//SOLUTION BEGIN
void pre() throws Exception{}
void solve(int TC) throws Exception{
int[] lpf = new int[MX];
for(int i = 2; i*i< MX; i++)
if(lpf[i]==0)
for(int j = i*i; j< MX; j+= i)
if(lpf[j]==0)
lpf[j] = i;
for(int i = 0; i< MX; i++)if(lpf[i]==0)lpf[i] = i;
int n = ni();
int[] a = new int[2+n];a[0] = 1;a[1+n] = 1;
int[][][] preM = new int[2+n][7][], sufM = new int[2+n][7][];
preM[0] = new int[][]{};
sufM[n+1] = new int[][]{};
for(int i = 1; i<= n; i++)a[i] = ni();
long[] pre = new long[2+n], suf = new long[2+n];
for(int i = 1; i<= n+1; i++){
pre[i] = pre[i-1];
int mx = 0, cur = 1;
int pp = 0, cp = 0;
while(a[i]>1){
int x = lpf[a[i]];
while(a[i]%x==0)a[i]/=x;
cur*=x;
while(pp<preM[i-1].length && preM[i-1][pp][0]<x)pp++;
int y = 1;
if(pp<preM[i-1].length && preM[i-1][pp][0]==x)y+=preM[i-1][pp][1];
mx = Math.max(mx, y);
preM[i][cp++] = new int[]{x, y};
}
preM[i] = Arrays.copyOfRange(preM[i], 0, cp);
a[i] = cur;
pre[i] += mx;
}
for(int i = n; i>= 0; i--){
suf[i] = suf[i+1];
int mx = 0,cur = a[i];
int pp = 0, cp = 0;
while(cur>1){
int x = lpf[cur];
while(cur%x==0)cur/=x;
while(pp< sufM[i+1].length && sufM[i+1][pp][0]<x)pp++;
int y = 1;
if(pp<sufM[i+1].length && sufM[i+1][pp][0]==x)y+=sufM[i+1][pp][1];
mx = Math.max(mx, y);
sufM[i][cp++] = new int[]{x, y};
}
sufM[i] = Arrays.copyOfRange(sufM[i], 0, cp);
suf[i]+=mx;
}
long ans = 0;
long[][] tmp = new long[15][3];
long[] sum = new long[1<<12], product = new long[1<<12], prevMask = new long[1<<12];
int[] mp = new int[1<<12];int[] pow = new int[12];
for(int i = 0; i< 12; i++){mp[1<<i] = i;pow[i] = 1<<i;}
for(int i = 1; i< 1<<12; i++)mp[i] = Math.max(mp[i], mp[i-1]);
for(int i = 1; i<= n; i++){
int c = 0;
for(int j = 0, k = 0; j<preM[i-1].length || k < sufM[i+1].length; ){
int x,le,ri;
if(j< preM[i-1].length && (k==sufM[i+1].length || sufM[i+1][k][0] > preM[i-1][j][0])){
x = preM[i-1][j][0];le = i-preM[i-1][j][1]; ri = i;
j++;
}else if(k< sufM[i+1].length && (j==preM[i-1].length || preM[i-1][j][0] > sufM[i+1][k][0])){
x = sufM[i+1][k][0];le = i; ri = i+sufM[i+1][k][1];
k++;
}else{
x = preM[i-1][j][0];le = i-preM[i-1][j][1]; ri = i+sufM[i+1][k][1];
j++;k++;
}
tmp[c][0] = x;
tmp[c][1] = le;
tmp[c][2] = ri;
c++;
}
Arrays.sort(tmp, 0, c, (long[] i1, long[] i2) -> Long.compare(i1[1], i2[1]));
loop:for(int mask = 1; mask< 1<<c; mask++){
int mask2 = mask^pow[mp[mask]];
if(mask2==0){
int pos = mp[mask];
product[mask] = tmp[pos][0];
prevMask[mask] = tmp[pos][2];
sum[mask] = (i-tmp[pos][1]+1)*(tmp[pos][2]-(i-1));
}else{
sum[mask] = sum[mask2];
int pos = mp[mask];
product[mask] = product[mask2];
if(product[mask]<MX)product[mask] *= tmp[pos][0];
prevMask[mask] = prevMask[mask2];
if(prevMask[mask]<tmp[pos][2]){
sum[mask] += (i-tmp[pos][1]+1)*(tmp[pos][2]-prevMask[mask]);
prevMask[mask] = tmp[pos][2];
}
}
if(product[mask]<MX)ans = Math.max(ans, pre[i-1]+suf[i+1]+sum[mask]);
}
}
pn(ans);
}
//SOLUTION END
void hold(boolean b)throws Exception{if(!b)throw new Exception("Hold right there, Sparky!");}
long mod = (long)1e9+7, IINF = (long)1e18;
final int INF = (int)1e9, MX = (int)5e5+1;
DecimalFormat df = new DecimalFormat("0.00000000000");
double PI = 3.1415926535897932384626433832792884197169399375105820974944, eps = 1e-8;
static boolean multipleTC = false, memory = false;
FastReader in;PrintWriter out;
void run() throws Exception{
in = new FastReader();
out = new PrintWriter(System.out);
int T = (multipleTC)?ni():1;
//Solution Credits: Taranpreet Singh
pre();for(int t = 1; t<= T; t++)solve(t);
out.flush();
out.close();
}
public static void main(String[] args) throws Exception{
if(memory)new Thread(null, new Runnable() {public void run(){try{new GCDS().run();}catch(Exception e){e.printStackTrace();}}}, "1", 1 << 28).start();
else new GCDS().run();
}
long gcd(long a, long b){return (b==0)?a:gcd(b,a%b);}
int gcd(int a, int b){return (b==0)?a:gcd(b,a%b);}
int bit(long n){return (n==0)?0:(1+bit(n&(n-1)));}
void p(Object o){out.print(o);}
void pn(Object o){out.println(o);}
void pni(Object o){out.println(o);out.flush();}
String n()throws Exception{return in.next();}
String nln()throws Exception{return in.nextLine();}
int ni()throws Exception{return Integer.parseInt(in.next());}
long nl()throws Exception{return Long.parseLong(in.next());}
double nd()throws Exception{return Double.parseDouble(in.next());}
class FastReader{
BufferedReader br;
StringTokenizer st;
public FastReader(){
br = new BufferedReader(new InputStreamReader(System.in));
}
public FastReader(String s) throws Exception{
br = new BufferedReader(new FileReader(s));
}
String next() throws Exception{
while (st == null || !st.hasMoreElements()){
try{
st = new StringTokenizer(br.readLine());
}catch (IOException e){
throw new Exception(e.toString());
}
}
return st.nextToken();
}
String nextLine() throws Exception{
String str = "";
try{
str = br.readLine();
}catch (IOException e){
throw new Exception(e.toString());
}
return str;
}
}
}
Feel free to Share your approach, If it differs. Suggestions are always welcomed. 