Cod sursa(job #17368)

Utilizator airineivAirinei Vasile airineiv Data 15 februarie 2007 19:30:10
Problema Fractii Scor 30
Compilator cpp Status done
Runda Arhiva de probleme Marime 8.03 kb
#include "stdio.h"

bool isPrime(char *prim, int n)
{
	if(n<=2)
		return true;
	if(n%2==0)
		return false;
	int i = (n-1)/2;
	if ((prim[i >> 3] & (1 << (i & 7))) == 0)    
		return true;
	return false;
}
void ProcessDivisors(unsigned int r, int m, unsigned divizor[10], long long &k)
{
	unsigned int j=r-1;
	if((unsigned)m<divizor[0])
	{
		k=m;
		return;
	}
	while(divizor[j] >(unsigned) m)
		j--;
	r=j+1;
	switch(r)
	{
	case 1:
		{
			k = m - m/divizor[0];
			break;
		}
	case 2:
		{
			k = m - m/divizor[0] - m/divizor[1] + m / (divizor[0] * divizor[1]);
			break;
		}
	case 3:
		{
			int p1 = divizor[0], p2 = divizor[1], p3 = divizor[2];
			k = m - m/p1 - m/p2 - m/p3 + m/(p1*p2) + m/(p1*p3) + m/(p2*p3) - m/(p1*p2*p3);
			break;
		}
	case 4:
		{
			int p1 = divizor[0], p2 = divizor[1], p3 = divizor[2], p4 = divizor[3];
			k = m - m/p1 - m/p2 - m/p3 - m/p4 + m/(p1*p2) + m/(p1*p3) + m/(p1*p4) + m/(p2*p3) + m/(p2*p4) + m / (p3*p4) - m/(p1*p2*p3) - m/(p1*p2*p4) - m/(p1*p3*p4) - m/(p2*p3*p4) + m/(p1*p2*p3*p4);
			break;
		}
	case 5:
		{
			int p1 = divizor[0], p2 = divizor[1], p3 = divizor[2], p4 = divizor[3], p5 = divizor[4];
			k = m - m/p1 - m/p2 - m/p3 - m/p4 -m/p5 + m/(p1*p2) + m/(p1*p3) + m/(p1*p4)  + m/(p1*p5) + m/(p2*p3) + m/(p2*p4) + m/(p2*p5) + m / (p3*p4) + m/(p3*p5) + m/(p4*p5) - m/(p1*p2*p3) - m/(p1*p2*p4) - m/(p1*p2*p5) - m/(p1*p3*p4) - m/(p1*p3*p5) - m/(p1*p4*p5) - m/(p2*p3*p4) - m/(p2*p3*p5) - m/(p2*p4*p5) - m/(p3*p4*p5);	
			k = k + m/(p1*p2*p3*p4) + m/(p1*p2*p3*p5) + m/(p1*p2*p4*p5) + m/(p1*p3*p4*p5) + m/(p2*p3*p4*p5) - m/(p1*p2*p3*p4*p5);
			break;
		}
	case 6:
		{
			int p1 = divizor[0], p2 = divizor[1], p3 = divizor[2], p4 = divizor[3], p5 = divizor[4], p6 = divizor[5];
			k = m - m/p1 - m/p2 - m/p3 - m/p4 - m/p5 - m/p6;
			k = k + m/(p1*p2) + m/(p1*p3) + m/(p1*p4)  + m/(p1*p5)  + m/(p1*p6) + m/(p2*p3) + m/(p2*p4) + m/(p2*p5) + m/(p2*p6) + m/(p3*p4) + m/(p3*p5) + m/(p3*p6) + m/(p4*p5) + m/(p4*p6) + m/(p5*p6);
			k = k - m/(p1*p2*p3) - m/(p1*p2*p4) - m/(p1*p2*p5) - m/(p1*p2*p6) - m/(p1*p3*p4) - m/(p1*p3*p5) - m/(p1*p3*p6) - m/(p1*p4*p5) - m/(p1*p4*p6) - m/(p1*p5*p6) - m/(p2*p3*p4) - m/(p2*p3*p5) - m/(p2*p3*p6) - m/(p2*p4*p5) - m/(p2*p4*p6) - m/(p2*p5*p6) - m/(p3*p4*p5) - m/(p3*p4*p6) - m/(p3*p5*p6) - m/(p4*p5*p6); 
			k = k + m/(p1*p2*p3*p4) + m/(p1*p2*p3*p5) + m/(p1*p2*p3*p6) + m/(p1*p2*p4*p5) + m/(p1*p2*p4*p6) + m/(p1*p2*p5*p6) + m/(p1*p3*p4*p5) + m/(p1*p3*p4*p6) + m/(p1*p3*p5*p6) + m/(p1*p4*p5*p6) + m/(p2*p3*p4*p5) + m/(p2*p3*p4*p6) + m/(p2*p3*p5*p6) + m/(p2*p4*p5*p6) + m/(p3*p4*p5*p6);
			k = k - m/(p1*p2*p3*p4*p5) - m/(p1*p2*p3*p4*p6) - m/(p6*p2*p3*p4*p5) - m/(p1*p6*p3*p4*p5) - m/(p1*p2*p6*p4*p5) - m/(p1*p2*p3*p6*p5) + m/(p1*p2*p3*p4*p5*p6);
			break;
		}

	case 7:
		{
			const int C71[7] = {0, 1, 2, 3, 4, 5, 6};
			const int C72[21][2] = {
				{0, 1},
				{0, 2},
				{0, 3},
				{0, 4},
				{0, 5},
				{0, 6},
				{1, 2},
				{1, 3},
				{1, 4},
				{1, 5},
				{1, 6},
				{2, 3},
				{2, 4},
				{2, 5},
				{2, 6},
				{3, 4},
				{3, 5},
				{3, 6},
				{4, 5},
				{4, 6},
				{5, 6}
			};
			const int C73[35][3] = {
				{0, 1, 2},
				{0, 1, 3},
				{0, 1, 4},
				{0, 1, 5},
				{0, 1, 6},
				{0, 2, 3},
				{0, 2, 4},
				{0, 2, 5},
				{0, 2, 6},
				{0, 3, 4},
				{0, 3, 5},
				{0, 3, 6},
				{0, 4, 5},
				{0, 4, 6},
				{0, 5, 6},
				{1, 2, 3},
				{1, 2, 4},
				{1, 2, 5},
				{1, 2, 6},
				{1, 3, 4},
				{1, 3, 5},
				{1, 3, 6},
				{1, 4, 5},
				{1, 4, 6},
				{1, 5, 6},
				{2, 3, 4},
				{2, 3, 5},
				{2, 3, 6},
				{2, 4, 5},
				{2, 4, 6},
				{2, 5, 6},
				{3, 4, 5},
				{3, 4, 6},
				{3, 5, 6},
				{4, 5, 6}
			};
			const int C74[35][4] = {
				{0, 1, 2, 3},
				{0, 1, 2, 4},
				{0, 1, 2, 5},
				{0, 1, 2, 6},
				{0, 1, 3, 4},
				{0, 1, 3, 5},
				{0, 1, 3, 6},
				{0, 1, 4, 5},
				{0, 1, 4, 6},
				{0, 1, 5, 6},
				{0, 2, 3, 4},
				{0, 2, 3, 5},
				{0, 2, 3, 6},
				{0, 2, 4, 5},
				{0, 2, 4, 6},
				{0, 2, 5, 6},
				{0, 3, 4, 5},
				{0, 3, 4, 6},
				{0, 3, 5, 6},
				{0, 4, 5, 6},
				{1, 2, 3, 4},
				{1, 2, 3, 5},
				{1, 2, 3, 6},
				{1, 2, 4, 5},
				{1, 2, 4, 6},
				{1, 2, 5, 6},
				{1, 3, 4, 5},
				{1, 3, 4, 6},
				{1, 3, 5, 6},
				{1, 4, 5, 6},
				{2, 3, 4, 5},
				{2, 3, 4, 6},
				{2, 3, 5, 6},
				{2, 4, 5, 6},
				{3, 4, 5, 6}
			};
			const int C75[21][5] = {
				{0, 1, 2, 3, 4},
				{0, 1, 2, 3, 5},
				{0, 1, 2, 3, 6},
				{0, 1, 2, 4, 5},
				{0, 1, 2, 4, 6},
				{0, 1, 2, 5, 6},
				{0, 1, 3, 4, 5},
				{0, 1, 3, 4, 6},
				{0, 1, 3, 5, 6},
				{0, 1, 4, 5, 6},
				{0, 2, 3, 4, 5},
				{0, 2, 3, 4, 6},
				{0, 2, 3, 5, 6},
				{0, 2, 4, 5, 6},
				{0, 3, 4, 5, 6},
				{1, 2, 3, 4, 5},
				{1, 2, 3, 4, 6},
				{1, 2, 3, 5, 6},
				{1, 2, 4, 5, 6},
				{1, 3, 4, 5, 6},
				{2, 3, 4, 5, 6}
			};
			const int C76[7][6] = {
				{0, 1, 2, 3, 4, 5},
				{0, 1, 2, 3, 4, 6},
				{0, 1, 2, 3, 5, 6},
				{0, 1, 2, 4, 5, 6},
				{0, 1, 3, 4, 5, 6},
				{0, 2, 3, 4, 5, 6},
				{1, 2, 3, 4, 5, 6}
			};
			int i, val = 0;
			for(j=0; j<7; j++)
				val += m/divizor[C71[j]];
			k = m-val;
			val = 0;
			for(i=0; i<21; i++)
			{
				int temp = 1;
				for(j=0; j<2; j++)
					temp *= divizor[C72[i][j]];
				val += m/temp;
			}
			k += val;
			val = 0;
			for(i=0; i<35; i++)
			{
				int temp = 1;
				for(j=0; j<3; j++)
					temp *= divizor[C73[i][j]];
				val += m/temp;
			}
			k -= val;
			val = 0;
			for(i=0; i<35; i++)
			{
				int temp = 1;
				for(j=0; j<4; j++)
					temp *= divizor[C74[i][j]];
				val += m/temp;
			}
			k += val;
			val = 0;
			for(i=0; i<21; i++)
			{
				int temp = 1;
				for(j=0; j<5; j++)
					temp *= divizor[C75[i][j]];
				val += m/temp;
			}
			k -= val;
			val = 0;
			for(i=0; i<7; i++)
			{
				int temp = 1;
				for(j=0; j<6; j++)
					temp *= divizor[C76[i][j]];
				val += m/temp;
			}
			k += val;
			k = k - m/(divizor[0]*divizor[1]*divizor[2]*divizor[3]*divizor[4]*divizor[5]*divizor[6]);
			break;
		}
	default:
		{
			
			break;
		}
	}
}

int main(void)
{
	FILE *fin, *fout;
	if((fin = fopen("fractii.in", "r"))==NULL)
		return 0;
	if((fout = fopen("fractii.out", "w"))==NULL)
		return 0;
	int N, mod, r;
	unsigned long long n, fi;
	unsigned int divizor[10];
	fscanf(fin, "%d", &N);
	int *t = new int[N+1]; 
	n = N;
	

	const int MAXSIZE = N/2/8+1;
	char *prim = new char[MAXSIZE];
	for(int k=0; k<MAXSIZE; k++)
		prim[k] = 0;

	//p[i] == 0 if 2*i + 1 is prime

	int i, j;
	for (i = 1; ((i * i) << 1) + (i << 1) <= n; i += 1) 
	{      
		if ((prim[i >> 3] & (1 << (i & 7))) == 0) 
		{
			for (j = ((i * i) << 1) + (i << 1); (j << 1) + 1 <= n; j += (i << 1) + 1) 
			{
				prim[j >> 3] |= (1 << (j & 7));
			}
		}
	}
	
	for(i=0; i<=N; i++)
		t[i] = 0;
	t[1] = 1;
		
	for(i=2; i<=N; i++)
	{
		if(isPrime(prim, i))
		{
			t[i] = i-1;
			unsigned long long k = i;
			while(k<=N)
			{
				t[k] = (i-1)*((int)k/i);
				k *= i;
			}
		}
	}
	for(i=2; i<=N; i++)
	{
		if(t[i] == 0)
		{
			for(j=2; j<=i; j++)
			{
				if(isPrime(prim, j))
				{
					int p = i;
					r = 0;
					while(p%j==0)
					{
						p = p/j;
						r++;
					}
					if(r>0 && t[p] != 0 && p%j!=0)
					{
						t[i] = t[p]*(j-1);
						for(int k=1; k<=r-1; k++)
							t[i] *= j;
						break;
					}
				}
			}
		}
	}
	n=N;
	for(i=2; i<=N; i++)
	{
		fi = t[i];
		mod = N%i;
		if(mod == 0)
			fi *= N/i;
		else
		{
			fi *= N/i;
			if(isPrime(prim, i))
				fi += mod;
			else
			{
				if(mod==1)
					fi++;
				else
				{
					long long k=0;
					r=0;
					int p=i;	
					for(j=2; j<=i; j++)
					{
						if(p%j == 0 && isPrime(prim, j))
						{
							while(p%j==0)
								p = p/j;
							divizor[r++] = j;
							if(p==1)
								break;
						}
					}
					
					ProcessDivisors(r, mod, divizor, k);
					fi += k;
				}
			}
		}
		n+=fi;
	}

	fprintf(fout, "%llu", n);
	fclose(fin);
	fclose(fout);

	delete []t;
	t = NULL;
	delete []prim;
	prim = NULL;
	return 0;
}