SPOJ LCMSUM – LCM Sum

Problem

Problem Link – SPOJ LCMSUM

Given $n$, calculate the sum $LCM(1,n) + LCM(2,n) + … + LCM(n,n)$, where $LCM(i,n)$ denotes the Least Common Multiple of the integers $i$ and $n$.

Solution

I recently solved this problem and found the solution really interesting. You can find the formula that solves this problem on OEIS. I will show you the derivation here.

In order to solve this problem, you need to know about Euler Phi Function, finding Divisor using Sieve and some properties of LCM and GCD.

Define SUM

Let us define a variable SUM which we need to find.

$SUM = LCM(1,n) + LCM(2,n) + … + LCM(n,n)$ (take $LCM(n,n)$ to the other side for now )
$SUM – LCM(n,n) = LCM(1,n) + LCM(2,n) + … + LCM(n-1,n)$

We know that $LCM(n,n) = n$.

$SUM – n = LCM(1,n) + LCM(2,n) + … + LCM(n-1,n)$ ($eq 1$)

Reverse and Add

$SUM – n = LCM(1,n) + LCM(2,n) + … + LCM(n-1,n)$ (Reverse $eq 1$ to get $eq 2$)
$SUM – n = LCM(n-1,n) + LCM(n-2,n) + … + LCM(1,n)$ ( $eq 2$ )

No what will we get if we add $eq 1$ with $eq 2$? Well, we need to do more work to find that.

Sum of $LCM(a,n) + LCM(n-a,n)$

$x = LCM(a,n) + LCM(n-a,n)$
$x = \frac{an}{gcd(a,n)} + \frac{(n-a)n}{gcd(n-a,n)}$ ($eq 3$)

Arghh. Now we need to prove that $gcd(a,n)$ is equal to $gcd(n-a,n)$.

If $c$ divides $a$ and $b$, then, $c$ will divide $a+b$ and $a-b$. This is common property of division. So, if $g = gcd(a,n)$ divides $a$ and $n$, then it will also divide $n-a$. Hence, $gcd(a,n)=gcd(n-a,n)$.

So $eq 3$ becomes:

$x = \frac{an}{gcd(a,n)} + \frac{(n-a)n}{gcd(a,n)}$
$x = \frac{an + n^2 -an}{gcd(a,n)}$.
$x = \frac{n^2}{gcd(a,n)}$.

Now, we can continue adding $eq 1$ with $eq 2$.

Reverse and Add Continued

$SUM – n = LCM(1,n) + LCM(2,n) + … + LCM(n-1,n)$ ($eq 1$)
$SUM – n = LCM(n-1,n) + LCM(n-2,n) + … + LCM(1,n)$ ( $eq 2$ Add them )$2(SUM-n)= \frac{n^2}{gcd(1,n)} + \frac{n^2}{gcd(2,n)} + … \frac{n^2}{gcd(n-1,n)}$
$2(SUM-n ) = \sum_{i=1}^{n-1}\frac{n^2}{gcd(i,n)}$
$2(SUM-n ) = n\sum_{i=1}^{n-1}\frac{n}{gcd(i,n)}$ ( take $n$ common )

Group Similar GCD

What are the possible values of $g = gcd(i,n)$? Since $g$ must divide $n$, $g$ needs to be a divisor of $n$. So we can list the possible values of $gcd(i,n)$ by finding the divisors of $n$.

Let $Z = n\sum_{i=1}^{n-1}\frac{n}{gcd(i,n)}$. Now, every time $gcd(i,n) = d$, where $d$ is a divisior of $n$, $n \times \frac{n}{d}$ is added to $Z$. Therefore, we just need to find, for each divisor $d$, how many times $n \times \frac{n}{d}$ is added to $Z$.

How many values $i$ can we put such that $gcd(i,n) = d$? There are $\phi(\frac{n}{d})$ possible values of $i$ for which we get $gcd(i,n)=d$. Therefore:

$2(SUM-n ) = n\sum_{i=1}^{n-1}\frac{n}{gcd(i,n)}$
$2(SUM-n ) = n\sum_{d |n, d \neq n } \phi(\frac{n}{d}) \times \frac{n}{d}$ (but, $\frac{n}{d}$ is also a divisor of $n$ )
$2(SUM-n ) = n\sum_{d |n, d \neq 1 } \phi(d) \times d$ ( when $d = 1$, we get $\phi(1)times 1$ )
$2(SUM-n ) = n( \sum_{d |n } ( \phi(d) \times d ) -1 )$
$2(SUM-n ) = n \sum_{d |n } (\phi(d) \times d ) – n$
$2SUM – 2n + n = n\sum_{d |n } \phi(d) \times d$
$2SUM – n = n\sum_{d |n } \phi(d) \times d$

$$2SUM = n( \sum_{d |n } \phi(d) \times d )+ n \\
2SUM = n ( \sum_{d |n } ( \phi(d) \times d ) + 1 ) \\
\therefore SUM = \frac{n}{2} ( \sum_{d |n } ( \phi(d) \times d ) + 1 )$$

Using this formula we can solve the problem.

Code


#include <bits/stdc++.h>

#define FOR(i,x,y) for(vlong i = (x) ; i <= (y) ; ++i)

using namespace std;
typedef long long vlong;

vlong res[1000010];
vlong phi[1000010];

void precal( int n ) {
    // Calculate phi from 1 to n using sieve
    FOR(i,1,n) phi[i] = i;
    FOR(i,2,n) {
        if ( phi[i] == i ) {
        for ( int j = i; j <= n; j += i ) {
            phi[j] /= i;
            phi[j] *= i - 1;
        }
    }
 }

 // Calculate partial result using sieve
 // For each divisor d of n, add phi(d)*d to result array
 FOR(i,1,n){
    for ( int j = i; j <= n; j += i ) {
        res[j] += ( i * phi[i] );
    }
 }
}

int main () {
    precal( 1000000 );

    int kase;
    scanf ( "%d", &kase );

    while ( kase-- ) {
        vlong n;
        scanf ( "%lld", &n );

        // We already have partial result in res[n]
        vlong ans = res[n] + 1;
        ans *= n;
        ans /= 2;

        printf ( "%lld\n", ans );
    }

    return 0;
}

We need to precalculate partial result using a sieve for all values of $N$. Precalculation has a complexity of $O(N \ ln(N))$. After pre-calculation, for each $N$ we can answer in $O(1)$.

Reference

  1. OEIS - A051193

Modular Exponentiation

Problem

Given three positive integers $B, P$ and $M$, find the value of $B^P \mod M$.

For example, $B=2$, $P=5$ and $M=7$, then $B^P \mod M = 2^5 \mod 7 = 32 \mod 7 = 4$.

This problem is known as [ref1] Modular Exponentiation. In order to solve this problem, we need to have basic knowledge about [ref2] Modular Arithmetic.

Brute Force – Linear Solution $O(P)$

As usual, we first start off with a naive solution. We can calculate the value of $B^P \mod M$ in $O(P)$ time by running a loop from $1$ to $P$ and multiplying $B$ to result.

int bigmodLinear ( int b, int p, int m ) {
    int res = 1 % m;
    b = b % m;
    for ( int i = 1; i <= p; i++ ) {
        res = ( res * b ) % m;
    }
    return res;
}

A simple solution which will work as long as the value of $P$ is small enough. But for large values of $P$, for example, $P > 10^9$, this code will time out. Also note that in line $2$, I wrote $res = 1 \mod m$. Some people tend to write $res = 1$. This will produce a wrong result when the value of $P$ is $0$ and value of $M$ is $1$.

Divide and Conquer Approach - $O(log_2(P) )$

According to [ref3] Wiki,

A divide and conquer algorithm works by recursively breaking down a problem into two or more sub-problems of the same (or related) type (divide), until these become simple enough to be solved directly (conquer).

That is, instead of solving the big problem $B^P \mod M$, we will solve smaller problems of similar kind and merge them to get the answer to the big problem.

So how do we apply this D&C (Divide and Conquer) idea?

Suppose we are trying to find the value of $B^P$. Then three thing can happen.

  1. Value of P is 0 ( Base Case ): $B^0$ is $1$. This will be the base case of our recursion function. Once we reach this state, we no longer divide the problem into smaller parts.
  2. Value of P is Even: Since $P$ is even, we can say that $B^P = B^{\frac{P}{2}} \times B^{\frac{P}{2}}$. For example, $2^{32} = 2^{16} \times 2^{16}$, $3^6 = 3^3 \times 3^3$. Therefore, instead of calculating the value of $x = B^P$, if we find the value of $y = B^{\frac{P}{2}}$, then we can get the value of $x$ as $x = y \times y$.
  3. Value of P is Odd: In this case we can simply say that $B^P = B \times B^{P-1}$.

Using these three states, we are can formulate a D&C code.

int bigmodRecur ( int b, int p, int m ) {
    if ( p == 0 ) return 1%m; // Base Case

    if ( p % 2 == 0 ) { // p is even
        int y = bigmodRecur ( b, p / 2, m );
        return ( y * y ) % m; // b^p = y * y
    }
    else {
        // b^p = b * b^(p-1)
        return ( b * bigmodRecur ( b, p - 1, m ) ) % m;
    }
}

In line $2$, we have the base case. We return $1 \mod m$ in case value of $m$ is $1$. Next on line $4$ we check if $p$ is even or not. If it is even, we calculate $A^{\frac{P}{2}}$ and return after squaring it. In line $8$, we handle the case when $P$ is odd.

At each step we either divide $P$ by $2$ or subtract $1$ from $P$ and then divide it by $2$. Therefore, there can be at most $2 \times log_2(P)$ steps in D&C approach. Hence complexity is $O(log_2(P) )$.

A Better Solution

A better solution to this problem exists. By using the Repeated Squaring algorithm, we can find the Modular Exponentiation faster than D&C. Repeated Squaring Algorithm has the same complexity as D&C, but it is iterative, thus does not suffer from recursion overhead.

We need to know about bitwise operations before we learn Repeated Squaring algorithm. Since we did not cover this topic yet, I won't write about this approach now. Hopefully, once we cover bitwise operations I will write another post

Resource

  1. Wikipedia - Modular Exponentiation
  2. forthright48 - Introduction to Modular Arithmetic
  3. Wikipedia - Divide and conquer algorithms
  4. forthright48 - Introduction to Number Systems

Leading Digits of Factorial

Problem

Given an integer $N$, find the first $K$ leading digits of $N!$.

For example, for $N=10$ and $K=3$, then first $3$ leading digits of $10! = 3628800$ is $362$.

Finding leading digits uses concepts similar to [ref1] Number of Trailing Zeroes of Factorial.

Brute Force Solution

Finding the value of $N!$ and then printing the first $K$ digits is a simple but slow solution. Using long long we can calculate value $N!$ up to $N \leq 20$ and using Big Integer we can calculate arbitrary $N!$ but with complexity worse than $O(N^2)$.

Solution Using Logarithm

In [ref1], we say that a logarithm of value $x$ is $y$ such that $x = 10^y$. For now let us find out leading digits of a value $x$ instead of $N!$. We will extend it to cover factorials later.

So, we know that $log_{10}(x) = y$, where $y$ will be some fraction. Let us separate $y$ into its integer and decimal part and call them $p,q$. For example, if $y = 123.456$, then $p=123$ and $q=0.456$.

Therefore, we can say that $log_{10}(x) = p + q$. Which means, $x = 10^y = 10^{p+q}=10^p \times 10^q$.

Now expand the values of $10^p$ and $10^q$. If $A=10^p$, then $A$ will simply be a power of $10$ since $p$ is an integer. To be more exact, $A$ will be $1$ with $p$ trailing zeroes. For example, $A=10^3 = 1000$. What about $B=10^q$?

Since $q$ is a fraction which is $0 \leq q < 1$, value of $B$ will be between $10^0 \leq B < 10^1$, that is, $1 \leq B < 10$.

Okay, we got the value of $A$ and $B$, what now? We know that if we multiply $A$ and $B$ we will get $x$. But don’t multiply them just yet. Think for a bit what will happen when we multiply a decimal number with $10$. If it is integer, it will get a trailing zero, e.g, $3 \times 10 = 30$. But if it is a fraction, its decimal point will shift to right, e.g $23.65 \times 10 = 236.5$. Actually, decimal points shifts for integer numbers too, since integer numbers are real numbers with $0$ as fraction, e.g $3 = 3.00$. So in either case multiplying $10$ shifts decimal point to the right.

So what happens if we multiply, $A$, which is just $10^p$ to $B$? Since $A$ has $10$ in it $p$ times, the decimal point of $B$ will shift to right $p$ times. That is all $A$ does to $B$ is change its decimal point. It does not change the digits of $B$ in any way. Thus, $B$ contains all the leading digits of $x$.

For example,
$log_{10}(5420) = 3.7339993 = 3 + 0.7339993$.
$\therefore B = 10^0.7339993 = 5.4200$.

So, if we need first $K$ leading digits of $x$, we just need to multiply $B$ with $10^{K-1}$ and then throw away the fraction part. That is $res = \lfloor B \times 10^{K-1} \rfloor$. Why $10^{K-1}$ not just $10^K$? That’s because we already have $1$ leading digit present in $10^q$ before shifting it.

Extending to Factorial

It is easy to extend the idea above to $N!$. First we need to find out the value of $y=log_{10}(N!)$.

$y= log_{10}(N!)$
$y= log_{10}(N \times (N-1) \times (N-2) \times… \times 3 \times 2 \times 1 )$
$y= log_{10}(N) + log_{10}(N-1) + log_{10}(N-2) + … + log_{10}(2) + log_{10}(1)$

So we can simply find out the value of $y$ by running a loop from $1$ to $N$ and taking its log value.

After that we decompose $y$ into $p$, integer part and $q$, fraction part. The answer will be $\lfloor 10^q \times 10^{K-1}\rfloor$.

Code

const double eps = 1e-9;

// Find the first K digits of N!
int leadingDigitFact ( int n, int k ) {
    double fact = 0;

    // Find log(N!)
    for ( int i = 1; i <= n; i++ ) {
        fact += log10 ( i );
    }

    // Find the value of q
    double q = fact - floor ( fact+eps );

    double B = pow ( 10, q );

    // Shift decimal point k-1 \times
    for ( int i = 0; i < k - 1; i++ ) {
        B *= 10;
    }

    // Don't forget to floor it
    return floor(B+eps);
}

The code does exactly what we discussed before. But note the $eps$ that we added when flooring value in line $12$ and $22$. This due to precision error when dealing with real numbers in C++. For example, due to precision error sometimes a value which is supposed to be $1$, becomes $0.9999999999999$. The difference between these two values is very small, but if we floor them both, the first one becomes $1$ whereas the second one becomes $0$. So in order to avoid this error, when flooring a positive value we add a small number ( epsilon = $0.000000001$ ) to the number.

Summary

We need to execute the following steps to find the first $K$ leading digits of a number $x$ ( in our problem $x=N!$ ):

  1. Find the log value of the number whose leading digits we are seeking. $y=log_{10}(x)$.
  2. Decompose $y$ into two parts. Integer part $p$ and fraction part $q$.
  3. The answer is $\lfloor 10^q \times 10^{K-1} \rfloor$.

Resource

  1. forthright48 - Number of Trailing Zeroes of Factorial