#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include "projecteuler.h"

int partition(void **array, int l, int r, int (*cmp)(void *lv, void *rv));

int is_prime(long int num)
{
    long int i, limit;

    if(num <= 3)
    {
        /* If num is 2 or 3 then it's prime.*/
        return num == 2 || num == 3;
    }

    /* If num is divisible by 2 or 3 then it's not prime.*/
    if(num % 2 == 0 || num % 3 == 0)
    {
        return 0;
    }

    /* Any number can have only one prime factor greater than its
     * square root. If we reach the square root and we haven't found
     * any smaller prime factors, then the number is prime.*/
    limit = floor(sqrt(num));

    /* Every prime other than 2 and 3 is in the form 6k+1 or 6k-1.
     * If I check all those value no prime factors of the number
     * will be missed. If a factor is found, the number is not prime
     * and the function returns 0.*/
    for(i = 5; i <= limit; i += 6)
    {
        if(num % i == 0 || num % (i + 2) == 0)
        {
            return 0;
        }
    }

    /* If no factor is found up to the square root of num, num is prime.*/
    return 1;
}

int is_palindrome(int num, int base)
{
    int reverse = 0, tmp;

    tmp = num;

    /* Start with reverse=0, get the rightmost digit of the number using
     * modulo operation (num modulo base), add it to reverse. Remove the
     * rightmost digit from num dividing num by the base, shift the reverse left
     * multiplying by the base, repeat until all digits have been inserted
     * in reverse order.*/
    while(tmp > 0)
    {
        reverse *= base;
        reverse += tmp % base;
        tmp /= base;
    }

    /* If the reversed number is equal to the original one, then it's palindrome.*/
    if(num == reverse)
    {
        return 1;
    }

    return 0;
}

/* Same function, using GMP Library for long numbers.*/
int is_palindrome_mpz(mpz_t num, int base)
{
    mpz_t tmp, reverse, rem;

    mpz_inits(tmp, reverse, rem, NULL);
    mpz_set(tmp, num);
    mpz_set_ui(reverse, 0);

    while(mpz_cmp_ui(tmp, 0) > 0)
    {
        mpz_mul_ui(reverse, reverse, base);
        mpz_tdiv_qr_ui(tmp, rem, tmp, base);
        mpz_add(reverse, reverse, rem);
    }

    if(!mpz_cmp(num, reverse))
    {
        mpz_clears(reverse, rem, NULL);
        return 1;
    }

    mpz_clears(reverse, rem, NULL);

    return 0;
}

long int gcd(long int a, long int b)
{
    /* Euclid's algorithm for the greatest common divisor:
     * gcd(a, 0) = a
     * gcd(a, b) = gcd(b, a modulo b)*/
    if(b == 0)
    {
        return a;
    }

    return gcd(b, a%b);
}

/* Least common multiple algorithm using the greatest common divisor.*/
long int lcm(long int a, long int b)
{
    return a * b / gcd(a, b);
}

/* Recursive function to calculate the least common multiple of more than 2 numbers.*/
long int lcmm(long int *values, int n)
{
    int i;
    long int value;

    /* If there are only two numbers, use the lcm function to calculate the lcm.*/
    if(n == 2)
    {
        return lcm(values[0], values[1]);
    }
    else
    {
        value = values[0];

        for(i = 0; i < n - 1; i++)
        {
            values[i] = values[i+1];
        }

        /* Recursively calculate lcm(a, b, c, ..., n) = lcm(a, lcm(b, c, ..., n)).*/
        return lcm(value, lcmm(values, n-1));
    }
}

/* Function implementing the Sieve or Eratosthenes to generate
 * primes up to a certain number.*/
int *sieve(int n)
{
    int i, j, limit;
    int *primes;

    if((primes = (int *)malloc(n*sizeof(int))) == NULL)
    {
        return NULL;
    }

    /* 0 and 1 are not prime, 2 and 3 are prime.*/
    primes[0] = 0;
    primes[1] = 0;
    primes[2] = 1;
    primes[3] = 1;

    /* Cross out (set to 0) all even numbers and set the odd numbers to 1 (possible prime).*/
    for(i = 4; i < n - 1; i += 2)
    {
        primes[i] = 0;
        primes[i+1] = 1;
    }

    /* If i is prime, all multiples of i smaller than i*i have already been crossed out.
     * if i=sqrt(n), all multiples of i up to n (the target) have been crossed out. So
     * there is no need check i>sqrt(n).*/
    limit = floor(sqrt(n));

    for(i = 3; i <= limit; i += 2)
    {
        /* Find the next number not crossed out, which is prime.*/
        if(primes[i])
        {
            /* Cross out all multiples of i, starting with i*i because any smaller multiple
             * of i has a smaller prime factor and has already been crossed out. Also, since
             * i is odd, i*i+i is even and has already been crossed out, so multiples are
             * crossed out with steps of 2*i.*/
            for(j = i * i; j < n; j += 2 * i)
            {
                primes[j] = 0;
            }
        }
    }

    return primes;
}

int count_divisors(int n)
{
    int i, limit, count = 0;

    /* For every divisor below the square root of n, there is a corresponding one
     * above the square root, so it's sufficient to check up to the square root of n
     * and count every divisor twice. If n is a perfect square, the last divisor is
     * wrongly counted twice and must be corrected.*/
    limit = floor(sqrt(n));

    for(i = 1; i <= limit; i++)
    {
        if(n % i == 0)
        {
            count += 2;
        }
    }

    if(n == limit * limit)
    {
        count--;
    }

    return count;
}

int find_max_path(int **triang, int n)
{
    int i, j;

    /* Start from the second to last row and go up.*/
    for(i = n - 2; i >= 0; i--)
    {
        /* For each element in the row, check the two adjacent elements
         * in the row below and sum the larger one to it. At the end,
         * the element at the top will contain the value of the maximum path.*/
        for(j = 0; j <= i; j++)
        {
            if(triang[i+1][j] > triang[i+1][j+1])
                triang[i][j] += triang[i+1][j];
            else
                triang[i][j] += triang[i+1][j+1];
        }
    }

    return triang[0][0];
}

int sum_of_divisors(int n, int proper)
{
    int i, sum = 1, limit;

    /* For each divisor of n smaller than the square root of n,
     * there is another one larger than the square root. If i is
     * a divisor of n, so is n/i. Checking divisors i up to square
     * root of n and adding both i and n/i is sufficient to sum
     * all divisors.*/
    limit = floor(sqrt(n));

    for(i = 2; i <= limit; i++)
    {
        if(n % i == 0)
        {
            sum += i;
            sum += n / i;
        }
    }

    /* If n is a perfect square, i=limit is a divisor and
     * has to be counted only once.*/
    if(n == limit * limit)
    {
        sum -= limit;
    }

    if(!proper)
    {
        sum += n;
    }

    return sum;
}

void swap(void **array, int i, int j)
{
    void *tmp;

    tmp = array[i];
    array[i] = array[j];
    array[j] = tmp;
}

void insertion_sort(void **array, int l, int r, int (*cmp)(void *lv, void *rv))
{
    int i, j;
    void *tmp;

    /* After this cycle the smallest element will be in the first position of the array.*/
    for(i = r; i > l; i--)
    {
        if(cmp(array[i], array[i-1]) < 0)
        {
            swap(array, i, i-1);
        }
    }

    /* For each element in the array (starting from i=2), move it to the left until a
     * smaller element on its left is found.*/
    for(i = l + 2; i <= r; i++)
    {
        tmp = array[i];
        j = i;

        while(cmp(tmp, array[j-1]) < 0)
        {
            array[j] = array[j-1];
            j--;
        }

        array[j] = tmp;
    }
}

int partition(void **array, int l, int r, int (*cmp)(void *lv, void *rv))
{
    int i = l -1, j = r;
    void *pivot;

    /* Arbitrarily selecting the rightmost element as pivot.*/
    pivot = array[r];

    while(1)
    {
        /* From the left, loop until an element greater than the pivot is found.*/
        while(cmp(array[++i], pivot) < 0);
        /* From the right, loop until an element smaller than the pivot is found
         * or the beginning of the array is reached.*/
        while(cmp(array[--j], pivot) > 0)
        {
            if(j == l)
            {
                break;
            }
        }

        /* If j<=i, array[j], which is smaller than pivot, is already on the left
         * of array[i], which is larger than pivot, so they don't need to be swapped.*/
        if(j <= i)
        {
            break;
        }

        /* If j>i, swap array[i] and array[j].*/
        swap(array, i, j);
    }

    /* Swap array[i] with pivot. All elements on the left of pivot are smaller, all
     * the elements on the right are bigger, so pivot is in the correct position
     * in the sorted array.*/
    swap(array, i, r);

    return i;
}

void quick_sort(void **array, int l, int r, int (*cmp)(void *lv, void *rv))
{
    int i;

    /* If the array is small, it's better to just use a simple insertion_sort algorithm.*/
    if(r - l <= 20)
    {
        insertion_sort(array, l, r, cmp);
        return;
    }

    /* Partition the array and recursively run quick_sort on the two partitions.*/
    i = partition(array, l, r, cmp);
    quick_sort(array, l, i-1, cmp);
    quick_sort(array, i+1, r, cmp);
}

/* Implements SEPA (Simple, Efficient Permutation Algorithm)
 * to find the next permutation.*/
int next_permutation(void **perm, int n, int (*cmp)(void *a, void *b))
{
    int i, key;

    /* Starting from the right of the array, for each pair of values
     * if the left one is smaller than the right, that value is the key.*/
    for(i = n - 2; i >= 0; i--)
    {
        if(cmp(perm[i], perm[i+1]) < 0)
        {
            key = i;
            break;
        }
    }

    /* If no left value is smaller than its right value, the
     * array is in reverse order, i.e. it's the last permutation.*/
    if(i == -1)
    {
        return -1;
    }

    /* Find the smallest value on the right of the key which is bigger than the key itself,
     * considering that the values at the right of the key are in reverse order.*/
    for(i = key + 1; i < n && cmp(perm[i], perm[key]) > 0; i++);

    /* Swap the value found and the key.*/
    swap(perm, key, i-1);
    /* Sort the values at the right of the key. This is
     * the next permutation.*/
    insertion_sort(perm, key+1, n-1, cmp);

    return 0;
}

int is_pandigital(int value, int n)
{
    int *digits;
    int i, digit;

    if((digits = (int *)calloc(n+1, sizeof(int))) == NULL)
    {
        fprintf(stderr, "Error while allocating memory\n");
        exit(1);
    }

    for(i = 0; i < n && value > 0; i++)
    {
        digit = value % 10;
        if(digit > n)
        {
            return 0;
        }
        digits[digit]++;
        value /= 10;
    }

    if(i < n || value > 0)
    {
        free(digits);
        return 0;
    }

    if(digits[0] != 0)
    {
        free(digits);
        return 0;
    }

    for(i = 1; i <= n; i++)
    {
        if(digits[i] != 1)
        {
            free(digits);
            return 0;
        }
    }

    free(digits);

    return 1;
}

int is_pentagonal(long int n)
{
    double i;

    /* A number n is pentagonal if p=(sqrt(24n+1)+1)/6 is an integer.
     * In this case, n is the pth pentagonal number.*/
    i = (sqrt(24*n+1) + 1) / 6;

    if(i == (int)i)
    {
        return 1;
    }
    else
    {
        return 0;
    }
}

/* Function implementing the iterative algorithm taken from Wikipedia
 * to find the continued fraction for sqrt(S). The algorithm is as
 * follows:
 *
 * m_0=0
 * d_0=1
 * a_0=floor(sqrt(n))
 * m_(n+1)=d_n*a_n-m_n
 * d_(n+1)=(S-m_(n+1)^2)/d_n
 * a_(n+1)=floor((sqrt(S)+m_(n+1))/d_(n+1))=floor((a_0+m_(n+1))/d_(n+1))
 * if a_i=2*a_0, the algorithm ends.*/
int *build_sqrt_cont_fraction(int i, int *period, int l)
{
    int mn = 0, mn1, dn = 1, dn1, a0, an, an1, count = 0, j;
    int *fraction;

    if((fraction = (int *)malloc(l*sizeof(int))) == NULL)
    {
        return NULL;
    }

    j = 0;
    a0 = floor(sqrt(i));
    an = a0;
    fraction[j] = an;
    j++;

    do
    {
        mn1 = dn * an - mn;
        dn1 = (i - mn1 * mn1)/ dn;
        an1 = floor((a0+mn1)/dn1);
        mn = mn1;
        dn = dn1;
        an = an1;
        count++;
        fraction[j] = an;
        j++;
    } while(an != 2 * a0);

    fraction[j] = -1;

    *period = count;

    return fraction;
}

/* Function to solve the Diophantine equation in the form x^2-Dy^2=1
 * (Pell equation) using continued fractions.*/
int pell_eq(int d, mpz_t x)
{
    int found = 0, j, period;
    mpz_t n1, n2, n3, d1, d2, d3, sol, tmp;
    int *fraction;

    /* Find the continued fraction for sqrt(d).*/
    if((fraction = build_sqrt_cont_fraction(d, &period, 100)) == NULL)
    {
        return -1;
    }

    /* Calculate the first convergent of the continued fraction.*/
    mpz_init_set_ui(n1, 0);
    mpz_init_set_ui(n2, 1);
    mpz_init_set_ui(d1, 1);
    mpz_init_set_ui(d2, 0);
    mpz_inits(n3, d3, sol, tmp, NULL);

    j = 0;
    mpz_mul_ui(n3, n2, fraction[j]);
    mpz_add(n3, n3, n1);
    mpz_mul_ui(d3, d2, fraction[j]);
    mpz_add(d3, d3, d1);
    j++;

    /* Check if x=n, y=d solve the equation x^2-Dy^2=1.*/
    mpz_mul(sol, n3, n3);
    mpz_mul(tmp, d3, d3);
    mpz_mul_ui(tmp, tmp, d);
    mpz_sub(sol, sol, tmp);

    if(mpz_cmp_ui(sol, 1) == 0)
    {
        mpz_set(x, n3);
        found = 1;
        free(fraction);
        mpz_clears(n1, n2, n3, d1, d2, d3, sol, tmp, NULL);
    }

    /* Until a solution is found, calculate the next convergent
     * and check if x=n and y=d solve the equation.*/
    while(!found)
    {
        mpz_set(n1, n2);
        mpz_set(n2, n3);
        mpz_set(d1, d2);
        mpz_set(d2, d3);
        mpz_mul_ui(n3, n2, fraction[j]);
        mpz_add(n3, n3, n1);
        mpz_mul_ui(d3, d2, fraction[j]);
        mpz_add(d3, d3, d1);

        mpz_mul(sol, n3, n3);
        mpz_mul(tmp, d3, d3);
        mpz_mul_ui(tmp, tmp, d);
        mpz_sub(sol, sol, tmp);

        if(mpz_cmp_ui(sol, 1) == 0)
        {
            mpz_set(x, n3);
            found = 1;
            free(fraction);
            mpz_clears(n1, n2, n3, d1, d2, d3, sol, tmp, NULL);
        }

        j++;

        if(fraction[j] == -1)
        {
            j = 1;
        }
    }

    return 0;
}

/* Function to check if a number is semiprime. Parameters include
 * pointers to p and q to return the factors values and a list of
 * primes.*/
int is_semiprime(int n, int *p, int *q, int *primes)
{
    int i, limit;

    /* If n is prime, it's not semiprime.*/
    if(primes[n])
    {
        return 0;
    }

    /* Check if n is semiprime and one of the factors is 2.*/
    if(n % 2 == 0)
    {
        if(primes[n/2])
        {
            *p = 2;
            *q = n / 2;
            return 1;
        }
        else
        {
            return 0;
        }
    }
    /* Check if n is semiprime and one of the factors is 3.*/
    else if(n % 3 == 0)
    {
        if(primes[n/3])
        {
            *p = 3;
            *q = n / 3;
            return 1;
        }
        else
        {
            return 0;
        }
    }

    /*Any number can have only one prime factor greater than its
    square root, so we can stop checking at this point.*/
    limit = floor(sqrt(n));

    /* Every prime other than 2 and 3 is in the form 6k+1 or 6k-1.
     * If I check all those value no prime factors of the number
     * will be missed. For each of these possible primes, check if
     * they are prime, then if the number is semiprime with using
     * that factor.*/
    for(i = 5; i <= limit; i += 6)
    {
        if(primes[i] && n % i == 0)
        {
            if(primes[n/i])
            {
                *p = i;
                *q = n / i;
                return 1;
            }
            else
            {
                return 0;
            }
        }
        else if(primes[i+2] && n % (i + 2) == 0)
        {
            if(primes[n/(i+2)])
            {
                *p = i + 2;
                *q = n / (i + 2);
                return 1;
            }
            else
            {
                return 0;
            }
        }
    }

    return 0;
}

/* If n=pq is semiprime, phi(n)=(p-1)(q-1)=pq-p-q+1=n-(p+4)+1
 * if p!=q. If p=q (n is a square), phi(n)=n-p.*/
int phi_semiprime(int n, int p, int q)
{
    if(p == q)
    {
        return n - p;
    }
    else
    {
        return n - (p + q) + 1;
    }
}

/* Function to calculate phi(n) for any n. If n is prime, phi(n)=n-1,
 * is n is semiprime, use the above function. In any other case,
 * phi(n)=n*prod(1-1/p) for every distinct prime p that divides n.*/
int phi(int n, int *primes)
{
    int i, p, q, limit;
    double ph = (double)n;

    /* If n is prime, phi(n)=n-1*/
    if(primes[n])
    {
        return n - 1;
    }

    /* If n is semiprime, use above function.*/
    if(is_semiprime(n, &p, &q, primes))
    {
        return phi_semiprime(n, p, q);
    }

    /* If 2 is a factor of n, multiply the current ph (which now is n)
     * by 1-1/2, then divide all factors 2.*/
    if(n % 2 == 0)
    {
        ph *= (1.0 - 1.0 / 2.0);

        do
        {
            n /= 2;
        } while(n % 2 == 0);
    }

    /* If 3 is a factor of n, multiply the current ph by 1-1/3,
     * then divide all factors 3.*/
    if(n % 3 == 0)
    {
        ph *= (1.0 - 1.0 / 3.0);

        do
        {
            n /= 3;
        } while(n % 3 == 0);
    }

    /*Any number can have only one prime factor greater than its
     * square root, so we can stop checking at this point and deal
     * with the only factor larger than sqrt(n), if present, at the end.*/
    limit = floor(sqrt(n));

    /* Every prime other than 2 and 3 is in the form 6k+1 or 6k-1.
     * If I check all those value no prime factors of the number
     * will be missed. For each of these possible primes, check if
     * they are prime, then check if the number divides n, in which
     * case update the current ph.*/
    for(i = 5; i <= limit; i += 6)
    {
        if(primes[i])
        {
            if(n % i == 0)
            {
                ph *= (1.0 - 1.0 / i);

                do
                {
                    n /= i;
                } while(n % i == 0);
            }
        }
        if(primes[i+2])
        {
            if(n % (i + 2) == 0)
            {
                ph *= (1.0 - 1.0 /(i + 2));

                do
                {
                    n /= (i + 2);
                } while(n % (i + 2) == 0);
            }
        }
    }

    /* After dividing all prime factors smaller than sqrt(n), n is either 1
     * or is equal to the only prime factor greater than sqrt(n). In this
     * second case, we need to update ph with the last prime factor.*/
    if(n > 1)
    {
        ph *= (1.0 - 1.0 / n);
    }

    return (int)ph;
}

/* Function implementing the partition function.*/
long int partition_fn(int n, long int *partitions, int mod)
{
    int k, limit;
    long int res = 0;

    /* The partition function for negative numbers is 0 by definition.*/
    if(n < 0)
    {
        return 0;
    }

    /* The partition function for zero is 1 by definition.*/
    if(n == 0)
    {
        partitions[n] = 1;

        return 1;
    }

    /* If the partition for the current n has already been calculated, return the value.*/
    if(partitions[n] != 0)
    {
        return partitions[n];
    }

    k = -ceil((sqrt(24*n+1)-1)/6);
    limit = floor((sqrt(24*n+1)+1)/6);

    while(k <= limit)
    {
        if(k != 0)
        {
            res += pow(-1, k+1) * partition_fn(n-k*(3*k-1)/2, partitions, mod);
        }
        k++;
    }

    /* Give the result modulo mod, if mod=!-1, otherwise give the full result.*/
    if(mod != -1)
    {
        partitions[n] = res % mod;

        return res % mod;
    }
    else
    {
        partitions[n] = res;

        return res;
    }
}

/* Function implementing Dijkstra's algorithm for a mxn matrix, where the matrix represents
 * a graph in which each value can be connected to two, three or four adjacent values.
 * The parameters are the graph matrix, a matrix storing the distances, the size m and n of
 * the matrix, two flags, up and back, that determines if the path can also go up and/or back
 * instead of only down and forward, and the starting row.*/
int dijkstra(int **matrix, int **distances, int m, int n, int up, int back, int start)
{
    int i, j, min_i, min_j, min;
    int **visited;

    if((visited = (int **)malloc(m*sizeof(int *))) == NULL)
    {
        return -1;
    }

    for(i = 0; i < m; i++)
    {
        if((visited[i] = (int *)calloc(n, sizeof(int))) == NULL)
        {
            return -1;
        }
    }

    /* Set the current distances to the maximum value.*/
    for(i = 0; i < m; i++)
    {
        for(j = 0; j < n; j++)
        {
            distances[i][j] = INT_MAX;
        }
    }

    /* Set the distance of the starting node to its value.*/
    i = start;
    j = 0;
    distances[i][j] = matrix[i][j];

    do
    {
        /* Visit the first node, and update the distance of its 2, 3 or 4
         * adjacent nodes.*/
        visited[i][j] = 1;

        if(i < m - 1 && distances[i][j] + matrix[i+1][j] < distances[i+1][j])
        {
            distances[i+1][j] = distances[i][j] + matrix[i+1][j];
        }

        if(up)
        {
            if(i > 0 && distances[i][j] + matrix[i-1][j] < distances[i-1][j])
            {
                distances[i-1][j] = distances[i][j] + matrix[i-1][j];
            }
        }

        if(j < n -1 && distances[i][j] + matrix[i][j+1] < distances[i][j+1])
        {
            distances[i][j+1] = distances[i][j] + matrix[i][j+1];
        }

        if(back)
        {
            if(j > 0 && distances[i][j] + matrix[i][j-1] < distances[i][j-1])
            {
                distances[i][j-1] = distances[i][j] + matrix[i][j-1];

            }
        }

        min = INT_MAX;

        /* Find the non visited node with the current minimum distance.*/
        for(i = 0; i < m; i++)
        {
            for(j = 0; j < n; j++)
            {
                if(!visited[i][j] && distances[i][j] <= min)
                {
                    min = distances[i][j];
                    min_i = i;
                    min_j = j;
                }
            }
        }

        i = min_i;
        j = min_j;
        /* Repeat until all nodes have been visited.*/
    } while(i != m - 1 || j != n - 1);

    return 1;
}

/* Function that calculates the radical of n, i.e. the product
 * of its distinct prime factors.*/
int radical(int n, int *primes)
{
    int i, limit, rad = 1;

    /* There can be only one prime factor of n greater than sqrt(n),
     * if we check up to sqrt(n) and divide all the prime factors we find,
     * at the end we will have 1 or the only prime factor larger than sqrt(n).*/
    limit = floor(sqrt(n));

    /* Check if n is divisible by two, and divide all 2s factors. Since 2
     * is the only even prime, we can then loop only on odd numbers.*/
    if(n % 2 == 0)
    {
        rad *= 2;

        do
        {
            n /= 2;
        } while(n % 2 == 0);
    }

    /* For each prime i, check if it's a factor of n, then divide n by i
     * until n % i != 0.*/
    for(i = 3; i <= limit; i+=2)
    {
        if(primes[i] && n % i == 0)
        {
            rad *= i;

            do
            {
                n /= i;
            } while(n % i == 0);
        }

        /* If n is prime, all other prime factors have been found.*/
        if(n == 1 || primes[n])
        {
            rad *= n;
            break;
        }
    }

    return rad;
}