# Problem

I have written code to solve this challenge:

Given a formula: GCD(x,y) means the GCD (Greatest Common Divisor) of x and y.

For example: if N=2, then:

• f(2) = 1/GCD(1, 4) + 2/GCD(2, 4) + 3/GCD(3, 4) + 4/GCD(4,4)
• f(2) = 1/1 + 2/2 + 3/1 + 4/4
• f(2) = 1 + 1 + 3 + 1
• f(2) = 6

Given N, find f(N).

Output:

Value of f(N) modulo 1.000.000.007 (109 + 7)

Edit: Output should be in modulo 10^9 + 7

I found a simplified formula using Wolframalpha, by finding the first 10 elements: (4^(n + 1) + 8) / 12

# Code

#include <bits/stdc++.h>
#define MOD 1000000007
using namespace std;
typedef short i16;
typedef unsigned short u16;
typedef int i32;
typedef unsigned int u32;
typedef long int i64;
typedef unsigned long int u64;
typedef float f32;
typedef double f64;

u64 pawa(u64 a, u64 b){
u64 tot = 1;
for (u64 i = 0; i < b; i++) {
tot *= a;
tot = tot % MOD;
}
}

int main(){
ios_base::sync_with_stdio(false);
cin.tie(NULL);
u64 n;
cin>>n;
cout<<(pawa(4, n + 1) + 8) / 12<<'\n';
return 0;
}


This code scales poorly - time limit exceeded. The upper corner constraint is 10^9. How do I improve my algorithm?

• Think about it a bit, or google it, but exponentiation should only use log(Exponent) multiplications, not Exponent ones. – Deduplicator Apr 11 at 9:54
• What's that MOD constant there for? Is there a part of the challenge that requires that (if so, it should be in the description). – Toby Speight Apr 11 at 12:20
• @TobySpeight Sorry, I forgot to write it, it's part of the challenge. Its included in the output description. – Andra Apr 11 at 12:35
• @Deduplicator, this particular case (a fixed base of 4) gets a nice easy head start by bit shifting before we even need to think about multiplication... – Toby Speight Apr 11 at 13:03

It's great that you found a closed form to generate the results. Shame that it was by search rather than by reasoning. (i/GCD(i, 2^n) is simply i shifted right until it's odd, so we get the sequence (in binary) 1**,** 10, 11, 100, 101, 110, 111, 1100, ... You can see that every other value is divided by 2, every fourth value divided by 4 and so on; that allows you to arrive at the closed form fairly quickly.

I'll note that the numerator and denominator have a common factor of 4, so could be reduced to $$\\frac{4^n+2}3\$$.

Now to the code:

<bits/stdc++.h> is non-portable, and extremely wasteful. Include only the headers that are needed; in this case just <iostream>.

Instead of using a preprocessor #define for a constant integer, always prefer to use a C++ constant. Macros don't respect context.

Avoid using namespace std;. It doesn't even save you typing, once the unnecessary stdio operations are removed.

The typedefs are misleading: the name u64 suggests that unsigned long happens to be 64 bits on your particular target, but it also suggests that you've written brittle code that depends on that assumption (and all the other typedefs appear to be unused). For this code, I think we really should be using std::uint_fast64_t from <cstdint>.

Moving on to main(): as hinted earlier, std::ios_base::sync_with_stdio() and std::cin::tie() are pointless, given that we're not using the C standard streams, so should be omitted. Also, there's no checking that streaming from std::cin was successful; that's easily fixed.

Now we need to look at pawa(). Given that we only ever call it with a equal to 4, then we can make a special-purpose function rather than a general modpow() - let's call it exp4mod() (meaning exponentiation base 4, modulo 1000000007).

For values of b less than 15, the result is less than 1000000007, so that is easily computed as 1u << (2*b). Values of 15 or more can be reduced using the identity $$\4^{14a+b} = (4^{14})^a4^b\$$, and then we can use a standard binary exponentiation for the first part and a shift for the second.

All that said, it's unclear why the exponentiation must be done modulo 1000000007 - the description says that the result be reduced modulo such a number. So perhaps there's a logic error there, and we really need to be reducing modulo 3000000021 (having reduced the fraction so we divide by 3)?

Here's the version I ended up with, having made the recommended changes (and modifying main() not to require input):

#include <cstdint>
#include <iostream>

using u64 = std::uint_fast64_t;

static u64 modpow(u64 x, u64 y, u64 mod)
{
u64 result = 1;
while (y) {
if (y%2) {
result = result * x % mod;
}
x = x * x % mod;
y /= 2;
}

return result;
}

static u64 exp4mod(u64 x)
{
static auto const mod = 3000000021;
static auto const chunk_size = 16;

static const u64 one = 1;
static auto const residue = (one << 2*chunk_size) % mod; // 4^15%mod = 73741817

const auto a = x / chunk_size;
const auto b = x % chunk_size;
return (modpow(residue, a, mod) << (2 * b)) % mod;
}

int main()
{
for (u64 n = 0;  n < 35;  ++n)
std::cout << n << ": " << (exp4mod(n) + 2) / 3 << '\n';
}

• Thanks for the answer. I just have read the modular exponentiation and solve it before, but I really appreciate the review! – Andra Apr 16 at 12:45