22

Estou precisando calcular o valor de uma função polinomial de grau arbitrário em valores arbitrários do domínio.

Gostaria saber como fazê-lo de modo que não seja ingênuo ao fazer o cálculo com ponto flutuante, considerando que no polinômio todos os coeficientes e valores de avaliação sejam relativamente normais e bem comportados (nada de NaN, infty, números denormalizados nem próximos da denormalização - não precisa se preocupar com essas validações).

  • Não se pode assumir que o polinômio será esparso (muitos coeficientes iguais a 0) ou que será cheio (todos os coeficientes diferentes de 0). Espera-se que a solução seja capaz de lidar com ambas as situações;

    • Exemplo de polinômio esparso: f(x)=2x^123-x^30+5.366x
    • Exemplo de polinômio cheio: f(x)=x^5-4x^4+3x^3-2x^2+x-9
  • Se seu algoritmo obter um comportamento absurdamente eficiente e preciso dependendo dessas condições, vale ressaltar que talvez, em um algoritmo adaptativo, valha a pena tentar rodar seu código; isto é, não precisa necessariamente utilizar o mesmo algoritmo para resolver polinômios esparsos e cheios, desde que a adaptação do algoritmo seja justificada;

  • Em casos de algoritmo adaptativo, espera-se que se tenha um algoritmo para validar a condição e seja demonstrado que validar a condição custa menos do que o ganho obtido de performance (ganho de precisão é difícil mensurar, mas aceito sugestões de medidas);

Gostaria de saber qual o erro esperado pela estratégia adotada, também. Por exemplo, resultados distintos podem ser obtidos se analisar da esquerda para a direita o polinômio f(x) = x^30+x^2+3000x se comparado com o polinômio f(x) = 3000x+x^2+x^30, mesmo que matematicamente sejam idênticos, podendo se agravar para valores muito grande de x (x >> 0). Como o seu algoritmo se comportará nesses casos?

Também gostaria de saber, levando em consideração a localidade de referência, a performance de cálculo desse algoritmo.

Casos de testes

  • Polinômio esparso: f(x)=2x^123-x^30+5.366x

    • f(x=1) = 6.366
    • f(x=10) = 2e123
    • f(x=0.01) = 0.05366
  • Polinômio cheio: f(x)=x^5-4x^4+3x^3-2x^2+x-9

    • f(x=1) = -10
    • f(x=10) = 62,801
    • f(x=0.01) = -8.9902

Não obtendo o valor exato, o erro deverá estar justificado e descrito se era previsto devido ao comportamento do algoritmo.

19
  • Relações de Girard para um polinômio de grau n. Commented 1/02/2019 às 5:06
  • 4
    @AugustoVasques, só que não me interessam raízes. Só conheço as relações de Girard aplicadas às raízes de polinômios. Commented 1/02/2019 às 5:08
  • 2
    Quem negativou não entendeu a sua pergunta. Só para esclarecer, a pergunta e sobre polinômio mínimo, o cara tem zilhões de entradas diárias em um website e para analisar essas dados não existe tempo de processamento disponível então ele coloca esses dados como uma matriz de coeficientes polinomiais e encontra o auto-vetor que é uma representatividade calculável dos dados originais, porém consome muito tempo de processamento. Commented 1/02/2019 às 12:17
  • 6
    Relações de Girard determinam as raízes do polinômio, mas e se 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...
    – Woss
    Commented 1/02/2019 às 12:33
  • 1
    Acho que fui rude. Porque relações de girard, porque tem um teorema de análise que diz que entre duas raízes reais consecutivas sempre existirá um ponto de inflexão e as derivadas da curva nessas estas raízes determina a magnitude da inflexão. Dotado dessas informações você pode aproximar a curva com um polinômio de Berstain. Commented 1/02/2019 às 13:41

2 Respostas 2

13

Resposta ainda incompleta, mas já tem conteúdo sobre polinômios completos.

Polinômios Completos

Por polinômios completos entendemos os polinômios que podem ser representados por

inserir a descrição da imagem aqui

Sendo ai diferente de zero; isto é, nenhum coeficiente é nulo e, portanto, o polinômio possuirá todos os termos relativos às potências de x.

Considerações Matemáticas

Tomemos como polinômio de exemplo um polinômio completo representado por

inserir a descrição da imagem aqui

Onde p(x) é o polinômio, ai é o coeficiente real não nulo de grau i e $x^i$ a i-ésima potência de x. Se expandirmos esta somatória, obteremos:

inserir a descrição da imagem aqui

Podemos observar que, excetuando-se o termo $0, todos os termos possuem x e sabemos que dada a propriedade associativa podemos colocar o termo em comum em evidência, obtendo:

inserir a descrição da imagem aqui

Se fizermos isso recursivamente com os termos dentro dos parenteses obteremos:

inserir a descrição da imagem aqui

Desta forma, para calcular o valor de p(x), podemos resolver:

  • bn = an
  • bn-1 = an-1 + bn x
  • bn-2 = an-2 + bn-1 x
  • ...
  • b1 = a1 + b2 x
  • b0 = a0 + b1 x

Resultando, assim, em p(x) = b0, de modo que reduzimos a resolução de um polinômio de grau n para a n resoluções de polinômios de primeiro grau. Esse é o método de Horner.

Comportamento

Enquanto para o cálculo de um ponto do polinômio de grau n no método tradicional demande n(n+1)/2 multiplicações e n adições, com o método de Horner o número de operações é reduzido para n multiplicações e n adições e sabemos que quanto menos operações não só o cálculo tende a ficar mais rápido como ficará menos suscetível à erros.

Uma implementação possível do método de Horner é, calculando p(x=1) = x5 - 4x4 + 3x3 - 2x2 + x - 9:

float x = 1;
float coefficients[] = {1, -4, 3, -2, 1, -9};
float result = 0;

for (int i = 0; i < sizeof(coefficients)/sizeof(float); i++) {
    result = coefficients[i] + result*x;
}

printf("%f\n", result);

Embora tenhamos o ganho com a redução de operações necessárias para o cálculo, ainda temos considerações a se fazer.

Primeiro, o problema intrínseco da computação em trabalhar com número com ponto flutuante. O tipo float da linguagem C é representado conforme a especificação IEEE 754 utilizando 32 bits, sendo 1 bit de sinal, 8 bits para representar o expoente e 23 bits para representar a mantissa, o que nos dá uma precisão de apenas 7 casas decimais. Ainda, quanto maior for a mantissa, menos casas decimais teremos disponíveis para a parte decimal, então espera-se que se o valor crescer demasiadamente a parte decimal se perca e, portanto, nosso resultado sofrerá truncamentos.

Segundo, o nosso resultado será definido por coefficients[i] + result*x. Da matemática, sabemos que a + b tenderá a a se a for muito maior que b, a >> b. Por exemplo, 10000 + 1 é aproximadamente 10000. Isso implicará que se result*x for muito maior que coefficients[i] o valor deste último será desprezível, fazendo também o truncamento do valor.

Para efeitos de comparação, ao executar o trecho de código acima, para x = 1, obtemos o resultado -10, que é o esperado. Já, quando executado para x = 1000, obtemos o resultado 996003050160128, enquanto o esperado era 996002998000991.

1
  • Acho que vale a pena ressaltar aqui que, se 1 > |x| > 0, a tendência é que o erro seja mínimo. Ok, está mais para conjectura minha, mas não vejo perder precisão devido aos "números grandes" Commented 6/02/2019 às 0:12
10
+50

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, 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
inserir a descrição da imagem aqui
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.

12
  • 1
    Bem lembrado da exponenciação do logaritmo, eu tinha me esquecido totalmente disso Commented 6/02/2019 às 14:43
  • Nota sobre a estratégia de usar o exponencial do logaritmo: causa uma perda de precisão bem chata =/ f(x) = (1.0 + x)^32.0 - e^(ln(1.0 + x) * 32.0) é igual a f(x) = 0 para x > -1. Entretanto, você acaba por perder alguma precisão. Por exemplo, para um polinômio esparso de grau 32, a estratégia do exp(log(x)*coef) já começa a perder precisão quando x >= 2. O Python usa as funções da libmath usadas pelo C/C++, então dá para extrapolar este exemplo para C. Mas isso não retira os méritos de sua resposta =D Commented 7/02/2019 às 11:54
  • Jefferson, no caso de trabalhar com double, sim. Porém em caso de trabalho com float pode-se internamente calcular logaritmo como double e obter exponencial convertida para float, o que dependendo de como fizer dá quase exatidão e ainda é possível ser bastante eficiente. Em outras palavras, isso depende de muitos fatores.
    – RHER WOLF
    Commented 7/02/2019 às 17:29
  • Tem uma estratégia que faz a exponenciação por expoente inteiro em o(log2(exp)), sendo exp o expoente do monômio. É semelhante à sua estratégia, porém memoiza dinamicamente ou resultados no lugar de usar um previamente definidos. Gostei muito da representação do polinômio esparso, estou vendo se consigo tornar isso em um algoritmo adaptativo, haha: Vou me debruçar no código para dar uma olhada em quantas operações de ponto flutuante foram feitas. Commented 8/02/2019 às 0:53
  • Sim, é aquela que trabalha com o código binário do inteiro para calcular, né? Ela faz o mesmo que citei acima, enxerga o expoente inteiro e como e0*2°+e1*2¹+e2*2²+e3*2³+... e portanto x^e como x^(e0+e1*2+e2*4+e3*8+...)=( x^1 )^e0*( x^2 )^e1*( x^4 )^e2*( x^8 )^e3*..., correto? A diferença desse com o que implementei é que esse recalcula as potências "quadrado em cima de quadrado", já o meu pré-calcula até certo ponto e não recalcula. Cada um tem suas vantagens.
    – RHER WOLF
    Commented 9/02/2019 às 11:56

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .