Code:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cmath>
#define GERAR_N_PONTOS 100000 // número de pontos que deseja amostrar
// na interseccao dos 9 conjuntos
#define PI 3.1415926535897932384626433832795028841968
/////////////////////////////////////////////
// funções auxiliares do programa
/////////////////////////////////////////////
inline
double parte_fracionaria( double d )
{
return fabs( d - (int)d );
}
//
// representação de um ponto no R7
//
class Ponto
{
public:
double x[7];
Ponto()
{
}
Ponto( const Ponto& p )
{
for( int i = 0; i < 7; i++ )
x[i] = p.x[i];
}
Ponto( const Ponto* p )
{
for( int i = 0; i < 7; i++ )
x[i] = p->x[i];
}
Ponto( double lambda, const Ponto *p )
{
for( int i = 0; i < 7; i++ )
x[i] = lambda * p->x[i];
}
double norma() const
{
return sqrt( norma2() );
}
double norma2() const
{
double sum = 0.0;
for( int i = 0; i < 7; i++ )
{
sum += x[i] * x[i];
}
return sum;
}
double produto_interno( const Ponto *x, const Ponto *y ) const
{
double sum = 0.0;
for( int i = 0; i < 7; i++ )
{
sum += x->x[i] * y->x[i];
}
return sum;
}
double distancia( const Ponto *x, const Ponto *y ) const
{
double sum = 0.0;
for( int i = 0; i < 7; i++ )
{
sum += ( x->x[i] - y->x[i] ) * ( x->x[i] - y->x[i] );
}
return sqrt( sum );
}
};
//
// representação de um conjunto no R7
//
class Conjunto
{
public:
virtual ~Conjunto() {}
virtual bool pertence( const Ponto *p ) = 0;
virtual Ponto* gera_ponto() { return NULL; };
double gera_valor_entre( double a, double b, double inclusive = true )
{
double comprimento = b - a;
double p;
if( inclusive )
{
p = (double)rand() / (double) RAND_MAX;
}
else
{
// exclui o máximo e o zero
p = rand() / ( RAND_MAX + 1 ) + 1;
}
return ( p * comprimento ) + a;
}
};
//
// representação do hipercubo (-1,1)^7 contido no R7
//
class Hipercubo : public Conjunto
{
public:
Hipercubo() {};
bool pertence( const Ponto *p )
{
for( int i = 0; i < 7; i++ )
if( fabs( p->x[i] ) > 1 )
return false;
return true;
}
Ponto* gera_ponto()
{
Ponto *p = new Ponto;
double f;
for( int i = 0; i < 7; i++ )
{
f = gera_valor_entre( -1.0, 1.0, true );
p->x[i] = f;
}
return p;
}
};
//
// representação da bola aberta de raio 1 (na norma 2) do R7
//
class Bola : public Conjunto
{
public:
Bola()
{
cubo = new Hipercubo();
};
bool pertence( const Ponto *p )
{
double norma = 0.0;
for( int i = 0; i < 7; i++ )
norma += p->x[i]*p->x[i];
return norma <= 1.0;
}
Ponto* gera_ponto()
{
Ponto *p = cubo->gera_ponto();
while( !pertence( p ) )
{
delete p;
p = gera_ponto();
}
return p;
}
private:
Hipercubo *cubo;
};
/////////////////////////////////////////////
// conjuntos do EP
/////////////////////////////////////////////
//
// R1 (bola aberta)
//
class R1 : public Bola {};
//
// R2
//
class R2 : public Conjunto
{
public:
bool pertence( const Ponto *p )
{
double soma = 0.0;
for( int i = 0; i < 6; i++ )
soma += pow( fabs( p->x[i] ), (double) (i + 1) );
return p->x[6] >= ( -0.5 + soma );
}
};
//
// R3
//
class R3 : public Conjunto
{
public:
bool pertence( const Ponto *p )
{
double soma = 0.0;
for( int i = 0; i < 7; i++ )
soma += p->x[i];
return soma <= 1.0;
}
};
//
// R4
//
class R4 : public Conjunto
{
public:
bool pertence( const Ponto *p )
{
double soma = 0.0;
for( int i = 0; i < 6; i++ )
soma += pow( fabs( p->x[i] ), 1.0 + fabs( p->x[i] ) );
return p->x[6] <= ( 1.0 - ( soma / 10.0 ) );
}
};
//
// R5
//
class R5 : public Conjunto
{
public:
bool pertence( const Ponto *p )
{
return ( p->x[0] * p->x[1] < 0.49 );
}
};
//
// R6
//
class R6 : public Conjunto
{
public:
bool pertence( const Ponto *p )
{
double soma = 0.0;
for( int i = 0; i < 7; i++ )
soma += p->x[i] * p->x[i];
return soma > 0.1;
}
};
//
// R7
//
class R7 : public Conjunto
{
public:
bool pertence( const Ponto *p )
{
Ponto *y = gera_ponto_no_cone( p );
Ponto *z = projecao_linear( p, y );
double d = p->distancia( p, z );
delete y;
delete z;
return d > 0.001;
}
Ponto* gera_ponto_no_cone( const Ponto *p )
{
Ponto *y = new Ponto();
y->x[6] = 0.0;
for( int i = 0; i < 6; i++ )
{
y->x[i] = p->x[i];
y->x[6] += p->x[i] * p->x[i];
}
y->x[6] = sqrt( y->x[6] );
return y;
}
Ponto* projecao_linear( const Ponto *x, const Ponto *y )
{
// calcula a projeção linear de x sobre y
double lambda = x->produto_interno( x, y ) / y->norma2();
Ponto *p = new Ponto( lambda, x );
return p;
}
};
//
// R8
//
class R8 : public Conjunto
{
public:
bool pertence( const Ponto *p )
{
// testaremos componente a componente, até que uma falhe (ou não)
double a;
double b;
double h = 0.01;
// componente 1
{
// sempre há um ponto y na curva gama cuja componente y1 é tão próxima
// quanto se queira da componente x1 de qualquer x na bola unitária,
// pois x1 pertence a (-1,1) e y1 = sen( t ) pertence a [-1, 1]
// o que faço aqui é determinar o intervalo (a, b) ao qual t deve pertencer
// para que esta distância não ultrapasse 0.01.
// este t limita as demais componentes.
double ta;
double tb;
ta = p->x[0] - h;
tb = p->x[0] + h;
if( ta < -1.0 ) ta = -1.0;
if( tb > 1.0 ) tb = 1.0;
// nota: pode ser que a < 0 e b > 0, tomar cuidado neste caso!
a = asin( ta );
b = asin( tb );
if( a > b )
{
printf( "R8::percente():: DEU MERDA -> a > b!!\n" );
exit( 1 );
}
}
// componentes 2 a 7
{
double c;
double d;
// a componente yi = sen( i*t + pi/7 )
for( int i = 1; i < 7; i++ )
{
// vou encontrar o intervalo (c, d) ao qual t deve pertencer para que
// a diferença das componentes seja inferior a 0.01.
double ta;
double tb;
double t1;
double t2;
ta = p->x[i] - h;
tb = p->x[i] + h;
if( ta < -1.0 ) ta = -1.0;
if( tb > 1.0 ) tb = 1.0;
t1 = asin( ta );
t2 = asin( tb );
c = ( t1 - PI/7.0 ) / (double) (i + 1);
d = ( t2 - PI/7.0 ) / (double) (i + 1);
if( c > d )
{
printf( "R8::percente():: DEU MERDA -> c > d!!\n" );
exit( 1 );
}
//printf( "[%i] -> (a, b) = (%2.6f, %2.6f) \t(c, d) = (%.6f, %.6f)\n", i, a, b, c, d );
// agora eu calculo a intersecção dos intervalos (c, d) e (a, b) e
// se for vazio, então há pelo menos uma componente xj, para j=1...i
// que difere de mais de 0.01 de yj, para qualquer y da curva gama.
if( c > a ) a = c;
if( d < b ) b = d;
if( b < a )
{
// a intersecção é vazia
return true;
}
}
}
return false;
}
};
//
// R9
//
class R9 : public Conjunto
{
public:
bool pertence( const Ponto *p )
{
double min = parte_fracionaria( 1000.0 * p->x[0] );
for( int i = 1; i < 7; i++ )
{
double atual = parte_fracionaria( 1000.0 * p->x[i] );
if( atual < min )
min = atual;
}
// o ponto não pode pertencer a este conjunto
return ! ( min < 0.01 );
}
};
//
// rotina para testar e analisar os conjuntos
//
void teste_e_analise()
{
unsigned int pontos_no_conjunto_R[9] = {0,0,0,0,0,0,0,0,0};
Conjunto* R[9] = { new R1(),
new R2(),
new R3(),
new R4(),
new R5(),
new R6(),
new R7(),
new R8(),
new R9()
};
Bola* bola = new Bola();
printf("\nGerando 10000 pontos na bola de raio 1...\n\n");
for( int i = 0; i < 10000; i++ )
{
Ponto *p = bola->gera_ponto();
bool temp1 = R[0]->pertence( p );
pontos_no_conjunto_R[0] += temp1;
bool temp2 = R[1]->pertence( p );
pontos_no_conjunto_R[1] += temp2;
bool temp3 = R[2]->pertence( p );
pontos_no_conjunto_R[2] += temp3;
bool temp4 = R[3]->pertence( p );
pontos_no_conjunto_R[3] += temp4;
bool temp5 = R[4]->pertence( p );
pontos_no_conjunto_R[4] += temp5;
bool temp6 = R[5]->pertence( p );
pontos_no_conjunto_R[5] += temp6;
bool temp7 = R[6]->pertence( p );
pontos_no_conjunto_R[6] += temp7;
bool temp8 = R[7]->pertence( p );
pontos_no_conjunto_R[7] += temp8;
bool temp9 = R[8]->pertence( p );
pontos_no_conjunto_R[8] += temp9;
delete p;
}
printf( "-- R1 possui %10i (rejeitou %10i) --\n", pontos_no_conjunto_R[0], 10000 - pontos_no_conjunto_R[0] );
printf( "-- R2 possui %10i (rejeitou %10i) --\n", pontos_no_conjunto_R[1], 10000 - pontos_no_conjunto_R[1] );
printf( "-- R3 possui %10i (rejeitou %10i) --\n", pontos_no_conjunto_R[2], 10000 - pontos_no_conjunto_R[2] );
printf( "-- R4 possui %10i (rejeitou %10i) --\n", pontos_no_conjunto_R[3], 10000 - pontos_no_conjunto_R[3] );
printf( "-- R5 possui %10i (rejeitou %10i) --\n", pontos_no_conjunto_R[4], 10000 - pontos_no_conjunto_R[4] );
printf( "-- R6 possui %10i (rejeitou %10i) --\n", pontos_no_conjunto_R[5], 10000 - pontos_no_conjunto_R[5] );
printf( "-- R7 possui %10i (rejeitou %10i) --\n", pontos_no_conjunto_R[6], 10000 - pontos_no_conjunto_R[6] );
printf( "-- R8 possui %10i (rejeitou %10i) --\n", pontos_no_conjunto_R[7], 10000 - pontos_no_conjunto_R[7] );
printf( "-- R9 possui %10i (rejeitou %10i) --\n\n", pontos_no_conjunto_R[8], 10000 - pontos_no_conjunto_R[8] );
for( int i = 0; i < 9; i++ )
delete R[i];
delete bola;
}
/////////////////////////////////////////////
// função principal do EP
/////////////////////////////////////////////
int main( int argc, char *argv[] )
{
// lista de conjuntos, do menor para o maior
//
// R2
// R3
// R9
// R7
// R6
// R4
// R1, R5, R8
//
// farei as comparações nesta ordem para minimizar os cálculos
// rodar esta rotina apenas para testar os conjuntos e amostrar as suas relevâncias
if( argc == 2 )
{
if( strcmp( argv[1], "testar" ) == 0 )
{
teste_e_analise();
return 1;
}
else if( strcmp( argv[1], "rodar" ) != 0 )
return 0;
}
else
{
printf( "\nForma de utilizacao do programa:\n" );
printf( " %s testar\t-- roda a rotina de teste dos conjuntos\n", argv[0] );
printf( " %s rodar\t-- calcula o volume da interseccao dos conjuntos\n\n", argv[0] );
return 0;
}
Conjunto* R[9] = { new R1(), new R2(), new R3(), new R4(), new R5(),
new R6(), new R7(), new R8(), new R9() };
Hipercubo* cubo = new Hipercubo();
int i, j;
printf( "\nGerando pontos do hipercubo (-1,1)^7 contido no R7\n ateh que %i destes pertencam aos 9 conjuntos...\n\n",
GERAR_N_PONTOS );
for( i = j = 0; i < GERAR_N_PONTOS; i++ )
{
bool ponto_dentro;
do
{
ponto_dentro = false;
Ponto *p = cubo->gera_ponto();
if( R[1]->pertence( p ) &&
R[2]->pertence( p ) &&
R[8]->pertence( p ) &&
R[6]->pertence( p ) &&
R[5]->pertence( p ) &&
R[3]->pertence( p ) &&
R[0]->pertence( p ) &&
R[4]->pertence( p ) &&
R[7]->pertence( p ) )
{
ponto_dentro = true;
}
if( j % 500000 == 0 )
{
printf( "-- %8i pontos gerados e %6i pertencentes a todos os conjuntos --\n", j, i );
}
delete p;
j++;
}
while( !ponto_dentro );
}
printf( "-- %8i pontos gerados e %6i pertencentes a todos os conjuntos --\n", j, i );
printf( "\n\n" );
double porcentagem = (double) i / (double) j * 100.0;
printf( "Pertencem aos 9 conjuntos: %i de %i [%.4f%%]\n", i, j, porcentagem );
printf( "Volume do hipercubo: 2^7 = 128\n" );
printf( "Volume da interseccao dos conjuntos: %.6f * 128 = %.4f\n\n",
porcentagem / 100.0, 128.0 * porcentagem / 100.0 );
for( int i = 0; i < 9; i++ )
delete R[i];
delete cubo;
return 1;
}