Polinômios Cheios
Em casos de polinômios cheios, creio ser mais recomendado estruturar o polinômio em forma de array de coeficientes, fazendo de cada célula um monômio componente do polinômio e para cada célula o índice corresponde ao expoente do parâmetro e o elemento lá armazenado, o coeficiente. Por exemplo, o array 5,4,-1,3
corresponde a (5)*x°+(4)*x¹+(-1)*x²+(3)*x³
, portanto 5+4x-x²+3x³
.
P(c0,c1,c2,...,cn-2,cn-1,cn)(x) = c0 + c1*x + c2*x² + ... + cn-2*xn-2 + cn-1*xn-1 + cn*xn
O motivo é que isso aproveita a característica de polinômio cheio em ao menos dois aspectos: (1) possibilita que seja necessário armazenar apenas um coeficiente por monômio e mais nenhum outro dado, nem sequer expoente e (2) permite calcular o valor do polinômio dado um argumento fazendo uso de apenas uma soma e uma multiplicação por monômio.
P(c0,c1,c2,...,cn-2,cn-1,cn)(x) =
cn*xn + cn-1*xn-1 + cn-2*xn-2 + ... + c2*x2 + c1*x1 + cn-0*x0 =
( cn*x + cn-1 )*xn-1 + cn-2*xn-2 + ... + c2*x2 + c1*x1 + cn-0*x0 =
P(c0,c1,c2,...,cn-2, cn*x + cn-1 )(x)
Percebe-se uma relação recursiva, né? Determinando um argumento para o parâmetro, pode-se com uma soma e uma multiplicação (cn*x+cn-1) reestruturar o polinômio de modo que se una os valores dos dois últimos monômios como se fossem um. Isso pode se repetir até que se reste apenas um valor de monômio, uma constante que é o resultado.
Polinômio( [c[0],c[1],...,c[n-1],c[n]] , x )
Se n<1 faça:
Se n<0 faça:
Retorna 0
Se n=0 faça:
Retorna c[0]
Retorna Polinômio( [c[0],c[1],...,c[n-1]+x*c[n]] , x )
Para que não faça chamada recursiva, pode-se implementar o seguinte algoritmo.
Polinômio( [c[0],c[1],...,c[n-1],c[n]] , x )
v <- 0
Para i=n,n-1,...,1,0 faça:
v <- v*x+c[i]
retorna v
Perceba que o segundo algoritmo somente tem uma soma e uma multiplicação num loop que executa n+1 vezes, ou seja, o número de monômios do polinômio. Sendo assim, o custo em número de operações matemáticas é O(n). Perceba também que em caso de coeficientes zerados na estrutura, mesmo representando um monômio nulo, é calculado também.
Polinômios Esparsos
Em casos de polinômios esparsos, acredito que a mesma solução adotada para polinômios cheios terá tanto o armazenamento desnecessariamente custoso quanto o desempenho desnecessariamente reduzido em casos específicos, afinal são muitos zeros que poderiam ser simplificados. Uma alternativa é estruturar com um array de monômios não nulos, onde cada monômio é armazenado com o expoente (e) da potência e o coeficiente (c). Por exemplo, (0,4),(9,3),(99,2),(999,1)
corresponde a 4+3x9+2x99+x999.
Pode-se fazer um algoritmo óbvio que calcula o valor de cada termo e os soma, sejam as potências calculadas usando produtos estratégicos ou função exponencial. Porém é possível, alternativamente, simular o segundo algoritmo de polinômios cheios de forma que se simplifique os excessivos produtos. Isso permite alcançar o resultado por meio de cálculos de potências com expoentes menores do que os dos monômios, o que requer menos produtos (mas se for via exponencial, dependendo de como ela é implementada, tanto faz). Além disso, se imediatamente associar um expoente a um meio de calcular a potência, pode ser que esse algoritmo seja adequado até em polinômios cheios.
P( (e0,c0) , (e1,c1) , (e2,c2) , ... , (en-2,cn-2) , (en-1,cn-1) , (en,cn) )(x) =
cn*xen + cn-1*xen-1 + cn-2*xen-2 + ... + c2*xe2 + c1*xe1 + cn-0*xe0 =
( cn*xen-en-1 + cn-1 )*xen-1 + cn-2*xen-2 + ... + c2*xe2 + c1*xe1 + cn-0*xe0 =
P( (e0,c0) , (e1,c1) , (e2,c2) , ... , (en-2,cn-2) , (en-1, cn*xen-en-1 + cn-1 ) )(x)
Em C/C++, é possível por meio de funções intrínsecas calcular logaritmos numa velocidade absurda e manualmente implementar exponenciais extremamente rápidas de tal maneira que x24 calculado pela fórmula exp(24*ln(x))
seja tão eficiente quanto x*=( x*=( x*=( x*=x*x ) ) )
ou até mais. De fato, a decisão da forma como se calcula as potências não é fácil. Há um universo de opções.
Algoritmo Escolhido
Decidiu-se utilizar uma abordagem baseada na simulação do segundo algoritmo voltado para polinômios cheios com potências realizadas via produtos. O algoritmo requer que o array de expoentes e coeficientes venha ordenado crescentemente por expoente (para que ei-ei-1 seja sempre não negativo) da mesma maneira que a versão para polinômios cheios (que a diferença é sempre +1).
Sabe-se que, refletindo escrita em base dois, um inteiro positivo pode ser retratado como um soma de potências de base dois e expoente inteiro (como 33=32+1
), portanto potências de expoente inteiro podem ser divididas em produtos de potências de mesma base e expoentes que são potência de dois (como x³³=x³²*x
). Tendo pré-calculada uma potência de expoente que é potência de dois, seu cálculo é economizado.
Vou assumir o fácil encaminhamento para o procedimento de cálculo de potência a partir de "switch-case". O algoritmo leva em conta que se começa sabendo que x°=1
, x¹=x
, x²=x*x
, x³=x*x²
e x^(4+n)=(x²)*(x²)*( x^n )
. Pode-se estender para o tanto que desejar (x^8
, x^16
, x^32
, etc), mas o código fica extenso e talvez complicado de entender, portanto vou restringir o conhecimento de x
, x²
e (x²)²
, bem como os procedimentos da maneira como descrevi.
Polinômio( [( e[0],c[0] ),( e[1],c[1] ),...,( e[n-1],c[n-1] ),( e[n],c[n] )] , x )
Se n<0 faça:
retorna 0
v <- c[n]
x2 <- x*x
x4 <- x2*x2
Para i=n-1,...,1,0 faça:
d <- e[i+1]-e[i]
switch d:
caso seja 0:
v <- v+c[i]
caso seja 1:
v <- v*x+c[i]
caso seja 2:
v <- v*x2+c[i]
caso seja 3:
v <- v*x2*x+c[i]
caso seja outro:
v <- v*x4
d <- d-4
vá para: switch d
d <- e[0]
switch d:
caso seja 0:
retorna v
caso seja 1:
retorna v*x
caso seja 2:
retorna v*x2
caso seja 3:
retorna v*x2*x
caso seja outro:
v <- v*x4
d <- d-4
vá para: switch d
retorna v
Para acelerar um pouco mais, pode-se fazer assim (conhecer cálculos definitivos não só até x^(4+i) mas até x^(8+i)).
Polinômio( [( e[0],c[0] ),( e[1],c[1] ),...,( e[n-1],c[n-1] ),( e[n],c[n] )] , x )
Se n<0 faça:
retorna 0
v <- c[n]
x2 <- x*x
x4 <- x2*x2
Para i=n-1,...,1,0 faça:
d <- e[i+1]-e[i]
switch d:
caso seja 0:
v <- v+c[i]
caso seja 1:
v <- v*x+c[i]
caso seja 2:
v <- v*x2+c[i]
caso seja 3:
v <- v*x2*x+c[i]
caso seja 4:
v <- v*x4+c[i]
caso seja 5:
v <- v*x4*x+c[i]
caso seja 6:
v <- v*x4*x2+c[i]
caso seja 7:
v <- v*x4*x2*x+c[i]
caso seja outro:
v <- v*x4*x4
d <- d-8
vá para: switch d
d <- e[0]
switch d:
caso seja 0:
retorna v
caso seja 1:
retorna v*x
caso seja 2:
retorna v*x2
caso seja 3:
retorna v*x2*x
caso seja 4:
retorna v*x4
caso seja 5:
retorna v*x4*x
caso seja 6:
retorna v*x4*x2
caso seja 7:
retorna v*x4*x2*x
caso seja outro:
v <- v*x4*x4
d <- d-8
vá para: switch d
retorna v
Perceba que, em relação à versão para polinômios cheios, ele adiciona dois produtos adicionalmente e utiliza um switch para selecionar um procedimento de cálculo que pode reduzir em até 75% (3/4) das operações de multiplicação. Se armazenar x^8, pode reduzir até 87,5% (7/8) e pode-se até mais.
Executando esse módulo em C++ que calcula os casos de teste fornecidos no enunciado usando esse algoritmo, tanto com floats quanto com doubles,
# include <iostream>
# include <string>
float p( float x , float *c , int *e , int n ){
float x2=x*x , x4=x2*x2 , v=c[n] ;
int i , d ;
for( i=n-1 ; i>=0 ; i-- ){
d = e[i+1]-e[i] ;
_1: switch( d ){
case 0: v = v+c[i] ; break ;
case 1: v = v*x+c[i] ; break ;
case 2: v = v*x2+c[i] ; break ;
case 3: v = v*x2*x+c[i] ; break ;
case 4: v = v*x4+c[i] ; break ;
case 5: v = v*x4*x+c[i] ; break ;
case 6: v = v*x4*x2+c[i] ; break ;
case 7: v = v*x4*x2*x+c[i] ; break ;
default: v = v*x4*x4 ; d-=8 ; goto _1 ;
}
}
d = e[0] ;
_2: switch( d ){
case 0: return v ;
case 1: return v*x ;
case 2: return v*x2 ;
case 3: return v*x2*x ;
case 4: return v*x4 ;
case 5: return v*x4*x ;
case 6: return v*x4*x2 ;
case 7: return v*x4*x2*x ;
default: v = v*x4*x4 ; d-=8 ; goto _2 ;
}
}
double p( double x , double *c , int *e , int n ){
double x2=x*x , x4=x2*x2 , v=c[n] ;
int i , d ;
for( i=n-1 ; i>=0 ; i-- ){
d = e[i+1]-e[i] ;
_1: switch( d ){
case 0: v = v+c[i] ; break ;
case 1: v = v*x+c[i] ; break ;
case 2: v = v*x2+c[i] ; break ;
case 3: v = v*x2*x+c[i] ; break ;
case 4: v = v*x4+c[i] ; break ;
case 5: v = v*x4*x+c[i] ; break ;
case 6: v = v*x4*x2+c[i] ; break ;
case 7: v = v*x4*x2*x+c[i] ; break ;
default: v = v*x4*x4 ; d-=8 ; goto _1 ;
}
}
d = e[0] ;
_2: switch( d ){
case 0: return v ;
case 1: return v*x ;
case 2: return v*x2 ;
case 3: return v*x2*x ;
case 4: return v*x4 ;
case 5: return v*x4*x ;
case 6: return v*x4*x2 ;
case 7: return v*x4*x2*x ;
default: v = v*x4*x4 ; d-=8 ; goto _2 ;
}
}
int main(){
float ch_cf[]={ -9 , 1 , -2 , 3 , -4 , 1 } , es_cf[]={ 5.366f , -1.0f , 2.0f } ;
double ch_cd[]={ -9 , 1 , -2 , 3 , -4 , 1 } , es_cd[]={ 5.366 , -1.0 , 2.0 } ;
int E[]={ 0 , 1 , 2 , 3 , 4 , 5 } , e[]={ 1 , 30 , 123 } ;
std::cout.precision(12) ;
std::cout << std::fixed << "PoliCh( 10.0 ) = " << p( 10 , ch_cf , E , 5 ) << " [" << p( 10 , ch_cd , E , 5 ) ;
std::cout << std::scientific << "] PoliEs( 10.0 ) = " << p( 10 , es_cf , e , 2 ) << " [" << p( 10 , es_cd , e , 2 ) << "]\n" ;
std::cout << std::fixed << "PoliCh( 1.00 ) = " << p( 1 , ch_cf , E , 5 ) << " [ " << p( 1 , ch_cd , E , 5 ) ;
std::cout << "] PoliEs( 1.00 ) = " << p( 1 , es_cf , e , 2 ) << " [" << p( 1 , es_cd , e , 2 ) << "]\n" ;
std::cout << "PoliCh( 0.01 ) = " << p( 0.01f , ch_cf , E , 5 ) << " [ " << p( 0.01 , ch_cd , E , 5 ) ;
std::cout << "] PoliEs( 0.01 ) = " << p( 0.01f , es_cf , e , 2 ) << " [" << p( 0.01 , es_cd , e , 2 ) << "]\n" ;
}
isto é impresso. Repare que em todos os casos se encontra um erro de aproximação relativo de menos de 0.000003% (3e-8, cerca de 2-25), o que é forte indicativo de que os erros são consequentes de arredondamentos de operações básicas, não erros de aproximação de fórmulas matemáticas.
PoliCh( 10.0 ) = 62801.000000000000 [62801.000000000000] PoliEs( 10.0 ) = inf [2.000000000000e+123]
PoliCh( 1.00 ) = -10.000000000000 [ -10.000000000000] PoliEs( 1.00 ) = 6.366000175476 [6.366000000000]
PoliCh( 0.01 ) = -8.990197181702 [ -8.990197039900] PoliEs( 0.01 ) = 0.053660001606 [0.053660000000]
É importante notar que não só PoliEs(10)
mas também PoliEs(0.01)
é caso de soma de números de magnitudes muito distintas. É como 1e99+1e9
, 1e99+(-1e9)
, -1e99+1e9
e -1e99+(-1e9)
que tem uma diferença de 90
de ordem de grandeza, portanto nem com vinte casas decimais se sente a modificação do número maior com a adição ou subtração do menor. Na prática, se arredonda o de menor valor absoluto para zero de tão desprezível em relação ao maior por não alcançar nem sequer o último dígito na precisão double.
Ainda assim, não se observam grandes erros nessas operações porque a mudança exata relativa do valor é ínfima ((( 1e99+1e9 )-( 1e99 ))/1e99 = 1e9/1e99 = 1e-90 ~ 0
), portanto desprezá-la proporciona erros insignificantes. Normalmente se observam erros é quando as magnitudes dos operados são mais equilibrados, de maneira que até possa afetar no último dígito armazenado.
Além disso, se observam erros maiores é quando a estrutura do valor é de natural baixa precisão (menos dígitos armazenados), os números envolvidos nas etapas do cálculo são "quebrados", há acúmulo de operações (portanto de arredondamentos, que geram erros) e erros matemáticos, como os erros de fórmulas que matematicamente aproximam outras, tipo séries de taylor abreviadas e interpolações.
Como o algoritmo utiliza todos os monômios fornecidos sem simplificações matemáticas e não usam fórmulas de aproximação (como potência usando aproximação exponencial), não se observam erros de motivação matemática. Os erros consequentes de limitações computacionais de armazenamento dos valores encontrados ao longo da execução podem ser observados em função do tipo de valor (float/double) usado e principalmente quando se utilizam polinômios cheios e grandes com coeficientes e argumento não representáveis com exatidão em binário.
Por exemplo, executando o código ajustado para o polinômio
encontramos, comparado aos valores exatos f(0.9)=1.80596532
e f(1.1)=-0.97304332
, com floats os valores f(0.9)=1.80596613884
e f(1.1)=-0.97304391861
(com, respectivamente, 8.2e-7
e 6e-7
de erros absolutos e 4.5e-7
e 6.2e-7
de erros relativos), o que indica erros em torno do 20º bit (entre 23). Já com doubles se observam os valores f(0.9)=1.80596531999999943
e f(1.1)=-0.97304332000000215
(com, respectivamente, 5.7e-16
e 2.2e-15
de erros absolutos e 3.2e-16
e 2.2e-15
de erros relativos), o que indica erro de f(0.9)
em torno do 51º bit (entre 52) e de f(1.1)
em torno do 48º bit.
Algoritmo Melhorado
Levando em conta ainda o comentário de Jefferson sobre o algoritmo o(log(e)) de potência, decidi reimplementar o algoritmo anterior usando-o em casos de expoentes grandes. Ele imprime o mesmo resultado.
# include <iostream>
# include <string>
float p( float x , float *c , int *e , int n ){
float x2=x*x , x4=x2*x2 ;
float xe[2]={1} , v=c[n] ;
int i , d , de ;
for( i=n-1 ; i>=0 ; i-- ){
d = e[i+1]-e[i] ;
_1: switch( d ){
case 0: v = v+c[i] ; break ;
case 1: v = v*x+c[i] ; break ;
case 2: v = v*x2+c[i] ; break ;
case 3: v = v*x2*x+c[i] ; break ;
case 4: v = v*x4+c[i] ; break ;
case 5: v = v*x4*x+c[i] ; break ;
case 6: v = v*x4*x2+c[i] ; break ;
case 7: v = v*x4*x2*x+c[i] ; break ;
default:
de = d>>3 ;
xe[1] = x4 ;
do {
xe[1] *= xe[1] ;
v *= xe[de&1] ;
} while( de>>=1 ) ;
d &= 7 ;
goto _1 ;
}
}
d = e[0] ;
_2: switch( d ){
case 0: return v ;
case 1: return v*x ;
case 2: return v*x2 ;
case 3: return v*x2*x ;
case 4: return v*x4 ;
case 5: return v*x4*x ;
case 6: return v*x4*x2 ;
case 7: return v*x4*x2*x ;
default:
de = d>>3 ;
xe[1] = x4 ;
do {
xe[1] *= xe[1] ;
v *= xe[de&1] ;
} while( de>>=1 ) ;
d &= 7 ;
goto _2 ;
}
}
double p( double x , double *c , int *e , int n ){
double x2=x*x , x4=x2*x2 ;
double xe[2]={1} , v=c[n] ;
int i , d , de ;
for( i=n-1 ; i>=0 ; i-- ){
d = e[i+1]-e[i] ;
_1: switch( d ){
case 0: v = v+c[i] ; break ;
case 1: v = v*x+c[i] ; break ;
case 2: v = v*x2+c[i] ; break ;
case 3: v = v*x2*x+c[i] ; break ;
case 4: v = v*x4+c[i] ; break ;
case 5: v = v*x4*x+c[i] ; break ;
case 6: v = v*x4*x2+c[i] ; break ;
case 7: v = v*x4*x2*x+c[i] ; break ;
default:
de = d>>3 ;
xe[1] = x4 ;
do {
xe[1] *= xe[1] ;
v *= xe[de&1] ;
} while( de>>=1 ) ;
d &= 7 ;
goto _1 ;
}
}
d = e[0] ;
_2: switch( d ){
case 0: return v ;
case 1: return v*x ;
case 2: return v*x2 ;
case 3: return v*x2*x ;
case 4: return v*x4 ;
case 5: return v*x4*x ;
case 6: return v*x4*x2 ;
case 7: return v*x4*x2*x ;
default:
de = d>>3 ;
xe[1] = x4 ;
do {
xe[1] *= xe[1] ;
v *= xe[de&1] ;
} while( de>>=1 ) ;
d &= 7 ;
goto _2 ;
}
}
int main(){
float ch_cf[]={ -9 , 1 , -2 , 3 , -4 , 1 } , es_cf[]={ 5.366f , -1.0f , 2.0f } ;
double ch_cd[]={ -9 , 1 , -2 , 3 , -4 , 1 } , es_cd[]={ 5.366 , -1.0 , 2.0 } ;
int E[]={ 0 , 1 , 2 , 3 , 4 , 5 } , e[]={ 1 , 30 , 123 } ;
std::cout.precision(12) ;
std::cout << std::fixed << "PoliCh( 10.0 ) = " << p( 10 , ch_cf , E , 5 ) << " [" << p( 10 , ch_cd , E , 5 ) ;
std::cout << std::scientific << "] PoliEs( 10.0 ) = " << p( 10 , es_cf , e , 2 ) << " [" << p( 10 , es_cd , e , 2 ) << "]\n" ;
std::cout << std::fixed << "PoliCh( 1.00 ) = " << p( 1 , ch_cf , E , 5 ) << " [ " << p( 1 , ch_cd , E , 5 ) ;
std::cout << "] PoliEs( 1.00 ) = " << p( 1 , es_cf , e , 2 ) << " [" << p( 1 , es_cd , e , 2 ) << "]\n" ;
std::cout << "PoliCh( 0.01 ) = " << p( 0.01f , ch_cf , E , 5 ) << " [ " << p( 0.01 , ch_cd , E , 5 ) ;
std::cout << "] PoliEs( 0.01 ) = " << p( 0.01f , es_cf , e , 2 ) << " [" << p( 0.01 , es_cd , e , 2 ) << "]\n" ;
}
Adicionalmente, para reduzir desperdícios em polinômios cheios pode-se adicionar ao algoritmo uma verificação do grau de esparsidade com operações adicionais de ordem O(1) a partir do maior grau de monômio e do número de termos para determinar qual é o procedimento mais adequado ao caso.
x
não for raiz? Ele quer analisar o polinômio em determinado ponto da curva, não verificar as raízes do mesmo. "[...] calcular o valor de uma função polinomial de grau arbitrário em valores arbitrários do domínio". A menos que a sua sugestão é utilizá-las para determinar as raízes, simplificar o polinômio e, consequentemente, simplificar o cálculo em valores arbitrários...