I tried solving the Coprime Triples problem. Got an AC at one go using the inclusion-exclusion principle. Then I learnt a little bit about the Mobius function and realized that this problem can be solved using it. I tried implementing my solution here and received a bunch of WAs. I really need some help figuring out what’s wrong with my code, given that it’s the first time I’m dealing with the Mobius function. Thank you. 
2 Likes
I tried making a very trivial test case generator (well documented) for this problem. It produces random vectors with random sizes.
I’m simply searching for the first mismatch of the values returned by the AC and WA code, so that I can work on an example (simple enough to solve manually) to figure out where exactly my WA code fails.
I guess my generator is too weak and too slow. It’d be really wonderful if anyone could help me improve the design or help me directly to find the error/s in the function MBSsolve.
Here's the code
/*~~~~~~~~~~~~~~~~~~~~~~~*/
/* */
/* */
/* */
/* ©raisinten */
/* Darshan Sen */
/* */
/* */
/* */
/*~~~~~~~~~~~~~~~~~~~~~~~*/
#include <bits/stdc++.h>
using namespace std;
#define endl "\n"
#ifdef _WIN32
#define gc getchar
#else
#define gc getchar_unlocked
#endif
#define ret return 0
#define what_is(x) cout << #x << " is: " << x << endl;
#define boost ios_base::sync_with_stdio(false), cin.tie(nullptr), cout.tie(nullptr)
#define MOD 1000000007
using ll = long long;
using llu = long long unsigned;
template <typename T>
void fs (T& x) {
x ^= x;
bool neg = false;
register char ch = gc();
while ((ch != '-') and (ch < '0' or ch > '9'))
ch = gc();
if (ch == '-') {
neg = true;
ch = gc();
}
for(; ch >= '0' and ch <= '9'; ch = gc())
x = (x << 1) + (x << 3) + ch - 48;
if (neg)
x *= -1;
}
template <typename T>
void showc (T b, T e) {
for (auto it = b; it != e; ++it)
cout << (*it) << " ";
cout << endl;
}
///////////////////////////////////////////////////////////
// max. bound for array elements
#define MAX_A 1'000'000
// seems to be working :)
int mobius[MAX_A + 1];
bool sieve[MAX_A + 1];
void calc_mobius (void) {
// init
fill(sieve, sieve + MAX_A + 1, true);
fill(mobius, mobius + MAX_A + 1, 1);
// mobius(non-square-free numbers) = 0
for (int i = 2, square; (square = i * i) <= MAX_A; ++i) {
for (int j = square; j <= MAX_A; j += square) {
mobius[j] = 0;
sieve[j] = false;
}
}
// mobius(the rest of the numbers) = (-1) ^ r
// r = number of prime factors
for (int i = 2; i <= MAX_A; ++i) {
if (sieve[i]) {
mobius[i] = -1;
for (int j = i << 1; j <= MAX_A; j += i) {
mobius[j] *= -1;
sieve[j] = false;
}
}
}
// mobius(1) = 1
mobius[1] = 1;
}
// AC solution below
int A[MAX_A + 1]; // stores the frequency of each element
ll res[MAX_A + 1]; // computes the number of triplets with gcd = index
int IEPsolve (const vector<int>& V, const int& N) {
// init
fill(A, A + MAX_A + 1, 0);
fill(res, res + MAX_A + 1, 0LL);
for (int _ = 0; _ < N; ++_)
++A[V[_]];
// inclusion
for (int i = 1; i <= MAX_A; ++i) {
ll temp = 0LLU;
for (int j = i; j <= MAX_A; j += i)
temp+= A[j];
res[i] = temp * (temp - 1) * (temp - 2) / 6;
}
// exclusion
for (int i = MAX_A; i >= 1; --i) {
for (int j = i << 1; j <= MAX_A; j += i) {
res[i] -= res[j];
}
}
// bye!
return res[1];
}
// WA solution :(
int frequency[MAX_A]; // stores frequency of index'th element
int multiples[MAX_A]; // stores number of multiples of the
// index'th element that exists in the
// input array
int MBSsolve (const vector<int>& V, const int& N) {
// init
fill(frequency, frequency + MAX_A, 0);
fill(multiples, multiples + MAX_A, 0);
for (int i = 0; i < N; ++i)
++frequency[V[i]];
// storing # of multiples
for (int i = 1; i <= MAX_A; ++i) {
for (int j = i; j <= MAX_A; j += i) {
multiples[i] += frequency[j];
}
}
// res = summation(mobius(i) x C(# of multiples of i, 3))
ll result = 0LL;
for (int i = 1; i <= MAX_A; ++i) {
ll n = multiples[i];
ll combinations = n * (n - 1) * (n - 2) / 6;
result += combinations * mobius[i];
}
// Ciao!
return result;
}
int main (void) {
#ifdef LOCAL
freopen("../test/input.txt", "r", stdin);
freopen("../test/output.txt", "w", stdout);
#endif
boost;
/*
Problem link:
https://www.codechef.com/problems/COPRIME3
Trying to solve COPRIME3 using Möbius function.
Gist: Find the number of triplets in an array
s.t. gcd of each is 1.
I've solved it using the inclusion-exclusion
principle. I extracted the main bits of the AC
code and placed it into the function IEPsolve.
Next I tried solving it using the Möbius function.
So, I made the function calc_mobius which fills
up the array mobius with the appropriate values
at respective indices using sieve. It has already been
tested against the code under the simple method here:
https://www.geeksforgeeks.org/program-mobius-function/
Seems to match 100%, which seems nice. :)
I then cut out the skeleton of my WA code,
which uses the Möbius function to compute
the number of triplets and put it into
the function named MBSsolve.
Below, I'm generating random vectors, named V
with size N and simply checking for the first
such vector, which causes a difference between
the results returned by the two functions.
I even tried running the program with an
infinite loop instead of placing a limit
to the number of iterations of the loop below.
My code simply isn't producing any result.
I'd really like to know how to improve this
test case generator.
*/
// seeding
srand(time(NULL));
// precalculating the mobius function
calc_mobius();
// limits for random test cases I generate below
int maxn = 10; // max. size of vector
int maxa = 100; // max. value of vector element
int iter = 20; // max. number of test cases
for (int _ = 0; _ < iter; ++_) {
// making random vector
int N = 1 + rand() % maxn;
vector<int> V(N);
for (int i = 0; i < N; ++i)
V[i] = 1 + rand() % maxa;
// storing results
int iep = IEPsolve(V, N);
int mbs = MBSsolve(V, N);
// check
if (iep != mbs) {
// gotcha!
cout << N << endl;
showc(V.begin(), V.end());
cout << iep << " " << mbs << endl;
break;
}
}
ret;
}
@ssjgz, I’ve noticed that you often create test case generators that are activated by the
--test flag. I’d really like to know how you would go about making one for this problem, given that both, the function we are testing and the one we are testing against are taking quite a bit of time to generate results. I know you’re a busy man, but it would be an honour to get your help, once you are free.
Thank you. 
1 Like
Not had a proper look at this yet, but this line might cause problems:
for (int i = 1; i <= MAX_A; ++i) {
ll n = multiples[i]; // <--
as multlples is declared as:
int multiples[MAX_A];
Edit:
Try changing:
int frequency[MAX_A]; // stores frequency of index'th element
int multiples[MAX_A]; // stores number of multiples of the
to
int frequency[MAX_A + 1]; // stores frequency of index'th element
int multiples[MAX_A + 1]; // stores number of multiples of the
and
fill(frequency, frequency + MAX_A, 0);
fill(multiples, multiples + MAX_A, 0);
to
fill(frequency, frequency + MAX_A + 1, 0);
fill(multiples, multiples + MAX_A + 1, 0);
This bug wouldn’t be caught by your test generator as it would require V to have an element with value MAX_A.
1 Like
facepalm How could I do that? Must be because of those sleepless nights XD.
Well done catching the error. Got an AC finally!
Thank you very much. 
But still, I’d really like to know about the testing part. What should be the approach for testing my solveOptimised quickly, where the solveBruteForce takes up a large amount of time? Thank you.
1 Like
This speeds up MBSsolve for small testcases quite a bit:
int MBSsolve (const vector<int>& V, const int& N) {
const int maxA = *max_element(V.begin(), V.end());
vector<int> frequency(maxA + 1, 0); // stores frequency of index'th element
vector<int> multiples(maxA + 1, 0); // stores number of multiples of the
// index'th element that exists in the
// input array
for (int i = 0; i < N; ++i)
++frequency[V[i]];
// storing # of multiples
for (int i = 1; i <= maxA; ++i) {
for (int j = i; j <= maxA; j += i) {
multiples[i] += frequency[j];
}
}
// res = summation(mobius(i) x C(# of multiples of i, 3))
ll result = 0LL;
for (int i = 1; i <= maxA; ++i) {
ll n = multiples[i];
ll combinations = n * (n - 1) * (n - 2) / 6;
result += combinations * mobius[i];
}
// Ciao!
return result;
}
(you can use arrays instead of the vectors, but I prefer vectors because a) dynamic arrays are not valid C++ (though you used static ones in yours); b) they always automatically initialise themselves; c) the -D_GLIBCXX_DEBUG compiler flag, with gcc, will detect out-of-bounds accesses etc which isn’t the case with arrays).
In general, I store the results of the solveBruteForce and just replay them into solveOptimised and check for discrepancies - then it matters a bit less how slow solveBruteForce is.
2 Likes
How did you get AC using Inclusion-Exclusion? If I am not wrong, it has Time Complexity O(2^{N}) for N sets. My idea is find out all the prime numbers that divide at-least one element of the input sequence, also for each of these primes find the number of numbers in the sequence that it divides.
We then apply the Inclusion-Exclusion principle to find the common multiples of two or more than two primes and subtract it from total numbers in the sequence. But this will make the complexity go exponential!
1 Like
Yes, that’s right, because none of the elements more than the max_element is going to have any contributions to the final result as the number of multiples they store is 0.
I actually came to know last month about the fact that variable length arrays don’t really work in C++, when I came across this article.
I guess I’ll shift to vectors permanently other than in some very rare situations, where, it is not required to initialize the array because static arrays are way faster in that case
(a tiny li’l subl instruction is all it takes :D).
I’m not sure I understand the usage of that flag. I wrote a program to test that out. It doesn’t seem to work any different from me excluding the flag.
Here's the program
#include <vector>
#define _GLIBCXX_DEBUG
int main (void) {
std::vector<int> V(3);
V[4] = 5;
return 0;
}
Result
It seems like it throws a Segmentation Fault only in the above case. I tried to change the index I’m trying to access @line 7, but to no avail.
How long does it take on an average to test your program completely on an average, before you get that AC? Do you check for the discrepancies using the diff command?
Thank you so much for the tips.
1 Like
Try moving the #define _GLIBCXX_DEBUG to the very top of the file - or even better, provide it on the command-line when compiling:
g++ -std=c++14 raisinten-test.cpp -g3 -O3 -D_GLIBCXX_DEBUG -Wall
Hmmm … not sure - I generally generate approx 10’000 small tests - depending on how quickly the brute force approach runs (some are O(N^2); some can be O(N!)), it generally doesn’t take more than a few mins to develop a testsuite.
I actually use a rather crap program of my own to do all the generation - I think others have posted equivalents from time-to-time, and they might be worth investigating.
It’s a bit crap because a) it only works under Linux + gcc and b) seems to have trouble feeding large testcases to the executable, which isn’t much of a problem in practice as I mainly use it for small ones.
General approach when starting a problem:
a) Copy over my Code Template; fill in the blanks; enable solveBruteForce; and implement the brute force solution.
b) Write a randomised testcase generator inside the body of the if (argc == 2 && string(argv[1]) == "--test") { ... }.
c) Run the testcase-generator.cpp with a regex to pipe in the generated testcases and extract the result of solveBruteForce - I actually have a small script which does most of this for me, including taking a copy of the .cpp file so I can work on it while it’s generating away 
d) Implement the optimised version, and use the --verify mode of testcase-generator.cpp to dump out the first testcase where it fails i.e. gives a different result to brute force.
e) When it passes all tests, I usually generate some large-ish random testcases (just as a kind of “smoke test” to make sure I haven’t done anything really stupid, performance-wise) and submit.
2 Likes
That’s an interesting approach you mentioned, but as you said, it’d receive a TLE.
My approach is quite similar.
After updating the array A with the frequencies of the respective elements, I store in the res array, the number of triplets that can be formed, s.t. the gcd of each such triplet is a multiple of the index. This is equal to ^{n}C_3, where, n is the number of multiples of the current index that are present in the array. If the size of the array is N, this can be done in O(N log N) by iterating over the multiples of an index like we do in a sieve.
Next comes the exclusion part. Remember, I said that I’m storing the number of triplets with a multiple of the index as the gcd. So, I clearly need to remove the contributions of the multiples of each index. This again, I perform in the same method discussed above (sieve) in O(N log N).
Now each index of res stores the number of triplets in the array with the gcd = index. Hence, I return the value stored at index 1 for the final result, i.e., the number of triplets with gcd = 1.
Hope this helps.
2 Likes
WTF! Finally used the flag correctly, thanks for the assistance. I’ve never seen anything like it before. I wonder what else the flag can do. Now I don’t need gdb or any IDE to debug my C++ code.
Wow, that’s a lot of work you’ve done. Thanks for sharing those files. Very interesting to look at. I’ll go through those more closely. Btw, great to know about your work process. Cheers!
1 Like