Problem 10: Square Free Function

Points:

Small Case: 5 Answer:139661

Large Case: 15 Answer: 7355940518480

A number is said to be square free if its prime decomposition contains no repeated factors. For example: 6=2*3 is square free, while 12=22 * 3 is not square free.
Find the sum of all N less than or equal to Z for which the value of f(N) is square-free. f(N) is defined as the smallest positive integer M such that AM = 1 (mod N) for all A co-prime to N.

e.g. If N=8
All numbers co-prime to 8 are 1, 3,5 and 7.
12 = 1 = 1 mod 8
32 = 9 = 1 mod 8
52 = 25 = 1 mod 8
72 = 49 = 1 mod 8
So 2 is the smallest such positive integer and since 2 is square-free, 8 will be included in our sum.

Small Case: Z = 1000
Large Case: Z = 10000000

Solution:

In number theory, the Carmichael function f(n) of a positive integer n, is defined as the smallest positive integer m such that am = 1 (mod n) for all a co-prime to n.

Mobius Function is defined by:
mobius(n) = 1 if n is a square-free positive integer with an even number of prime factors.
mobius(n) = −1 if n is a square-free positive integer with an odd number of prime factors.
mobius(n) = 0 if n is not square-free
and (Mobius(f(n)))2,the multiplicative function can be used to determine sqare free carmichael numbers.

Reference:Francesco Pappalardi,Filip Saidak and Igor E. Shparlinski,Square Free values of Carmichael Function, Journal of Number Theory 103(2003) 122–131

Code:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include<iostream>
using namespace std;
int A[10000000]={0},F[10000000]={0},isprime[3163];
int getF(int n, int p)
{
     int x,j,e,t;
     if (p==2)x=n/4;
     else x=(n/p)*(p-1);
     for (j=2;j*j<=x;j++)
     {
         if (isprime[j] && x%j==0)
         {
            e=0;
            t=x;
            while (t%j==0)
                  t/=j,e++;
            if (e>1)
               return 0;
         }
     }
     return 1;
}
int main()
{
    int d,n,e,p,k,i,j;
    for (n=2;n<3163;n++)
        isprime[n]=1;

    //Sieve for Eratosthenes for Prime
    //Storing the smallest prime which divides n.
    //If A[n]=0 it means it is prime number.
    for(d=2;d<3163;d++)
    {
        if(isprime[d])
        {
            for (n=d*d;n<3163;n+=d)
            {
                isprime[n]=0;
                A[n]=d;
            }
            for (;n<=10000000;n+=d)
                A[n]=d;
        }
    }

    //Applying the formula
    //F[n]=Mu^2(Lambda(n)) is Multiplicative
    F[1]=1;F[2]=1;F[3]=1;F[4]=1;
    for(n=5;n<=10000000;n++)
    {
        if (A[n]==0)
        {
           F[n]=getF(n,n);
        }
        else
        {
           p=A[n],k=n/p,e=p;
           while (k%p==0)
                 k/=p,e*=p;
           if (k==1) 
            F[n]=getF(n,p);
           else F[n]=F[k]*F[e];
        }
    }
    long long s=0;
    for (n=1;n<=10000000;n++)
        if(F[n])
                s+=n;
    cout<<s<<endl;
    return 0;
}