Procurando por um algoritmo arcsin

7

Alguém tem um algoritmo simples para calcular um arco-seno razoavelmente preciso? Por "simples", quero dizer algum tipo de polinômio que requer <= 5 multiplicado por amostra de saída. E por "razoavelmente preciso", quero dizer um algo cujo erro não exceda 10% quando o argumento de entrada estiver próximo de mais ou menos um. Pesquisei na web por um tempo, mas não encontrei nada imediatamente útil.

Richard Lyons
fonte
Isso pode dar algumas idéias stackoverflow.com/questions/5920467/...
geometrikal
Mas por que não apenas uma tabela de pesquisa?
geometrikal
Estou pensando em uma implementação em que a memória disponível é dolorosamente limitada. Portanto, não considerei nenhuma solução de 'tabela de consulta'. Obrigado por seus pensamentos.
Richard Lyons
Você permite raízes quadradas? Devido ao comportamento da função próxima a (inclinação infinita), uma aproximação polinomial não funciona bem lá. ±1
Yves Daoust
E o CORDIC, que leva apenas algumas adições e subtrações e sem multiplicações.
mattgately

Respostas:

5

Aqui está apenas uma versão polinomial :

arcsin(x)=x+12x33+1324x55+135246x77
function y = arcsin_test3(x)
    y = x.*(1+x.*x.*(1/6+ x.*x.*(3/(2*4*5) + x.*x.*((1*3*5)/(2*4*6*7)))))
endfunction

que parece ter cinco multiplicações (supondo que você possa salvar o resultado de x.*x) e três adições.

E o scilabenredo é:

insira a descrição da imagem aqui

Top é scilab's asinvs este, inferior é o erro entre os dois.


Resposta original

A raiz quadrada aqui pode ser um aborrecimento, mas pensei em escrevê-la porque parece divertida. :-)

Esta página sugere:

da página 81 do Handbook of Mathematics Functions, de Milton Abramowitz e Irene Stegun:

arcsin(x)=π/21x(a0+a1x+a2x2+a3x3),
que
a0=1.5707288a1=0.2121144a2=0.0742610a3=0.0187293

Eu implementei isso scilabe funciona bem, exceto em torno de . Apenas refletir sobre para uma aproximação muito melhor.x=10x11x0

O gráfico superior mostra scilaba asinfunção da comparação acima (em vermelho tracejado) em relação à minha mudança de verde.

O gráfico inferior mostra o erro da minha alteração (plotar isso e o original nos mesmos eixos significa que o verde parece zero em todos os lugares).

insira a descrição da imagem aqui

// 25770
function y = arcsin_test(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293

    xx = abs(x)

    y = %pi/2 - sqrt(1-x).*(a0 + a1*x + a2.*x.*x + a3.*x.*x.*x)

endfunction

function y = arcsin_test2(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293

    xx = abs(x)

    y = %pi/2 - sqrt(1-xx).*(a0 + a1*xx + a2.*xx.*xx + a3.*xx.*xx.*xx)

    y = y.*sign(x); 
endfunction

x = [-1: .0100001 : 1];

clf
subplot(211)
plot(x,arcsin_test2(x),'g.');
plot(x,arcsin_test(x),'r:');
plot(x,asin(x))
subplot(212)
//plot(x,(arcsin_test(x) - asin(x)),'r:')
plot(x,(arcsin_test2(x) - asin(x)),'g.')
Peter K.
fonte
2
"Handbook of Mathematical Functions" amor esse livro
geometrikal
11
Oh, dispara. Eu vi a página da web 'funções trigonométricas inversas' durante minha pesquisa na web, mas não rolei o suficiente para ver o material da 'série infinita'! Que vergonha. Peter K., esse é outro que devo a você. (O meu problema original para melhorar o desempenho de um diferencial digital de diferença central, que, creio eu, pode ser feito através da realização de uma operação arcsine.)
Richard Lyons
Sim, mas Rick, você não pode fazer uma série infinita. se for finito, os coeficientes ideais não serão exatamente o que você obtém ao truncar a série infinita. se você tem o MATLAB (eles têm uma licença de uso doméstico relativamente barata agora), posso enviar o código do MATLAB para a troca do Remez em função do desejo do seu coração.
22415 Robert bristow-johnson #
11
ATUALIZAÇÃO: Consegui resolver meu problema original (criando um diferenciador digital muito simples que melhorou o desempenho em comparação ao de um diferenciador de diferença central) sem usar a função arcsin (). Posso postar um blog em dsprelated.com descrevendo meus resultados. Agradeço a todos por sua ajuda!
Richard Lyons
3

Eu tenho uma implementação muito boa de aqui .arctan()

Eu acho que você pode usar a identidade:

arcsin(x)=arctan(x1x2)

para conseguir o que você quer.

Robert Bristow-Johnson
fonte
Seu link para várias funções é interessante, Robert. Apenas para rir, tentei implementar a função sqrt (1 + x). Em vez de usar os limites corretos de 0 a 4, estraguei tudo e usei 1 a 5. É claro que calculei um resultado incorreto. No entanto, quando eu multipliquei meu resultado "incorreto" por dois, acabei com o resultado correto. Interessante, hein?
Richard Lyons
Rick, tenho certeza de que as funções estão "corretas" (ou razoavelmente precisas) declaradas como estão com os limites de conforme indicado. para , é bom apenas para ; portanto, se você estiver entre e , precisará ter uma constante um pouco (a ) armazenada lá e você precisará saber a diferença entre um expoente par de e um expoente ímpar de . xx1x224222
Robert Bristow-Johnson
Eu também acredito que suas funções estão corretas. Eu estava apenas comentando um erro bobo de "limites de soma" que cometi e, ao cometer esse erro, calculei um resultado incorreto. Mas notei que meu resultado incorreto era exatamente a metade do resultado correto. Eu estava apenas dizendo que meu resultado incorreto tinha uma relação "interessante" com o resultado correto, só isso. Desculpe a confusão Robert.
Richard Lyons
2

A parte central da curva não é um problema real, pois é bastante linear e a aproximação de Taylor a dois ou três termos é um bom ponto de partida (o polinômio de mínimos quadrados se ajusta um pouco melhor).

Os lados são mais problemáticos por causa da inclinação infinita. Uma maneira de lidar é através da transformação

arcsin(x)=π2arcsin(1x2),

que envolve uma raiz quadrada.


Se seu argumento é representado com ponto flutuante, é obtida uma aproximação rápida da raiz quadrada pela metade do expoente e aplicando uma transformação linear à mantissa.z

Seja , com , então . Você pode aproximar por .z=m2e1m<2z=m2e/2m(21)(m+2)

  • desmonte o expoente (limpe-o produz a representação de );em
  • se for par, calcule ;e(21)(m+2)
  • se for ímpar, calcule ;e2(21)(m+2)
  • defina o expoente .e/2

insira a descrição da imagem aqui

Yves Daoust
fonte
Estou ansioso para experimentar esse interessante algoritmo de raiz quadrada!
Richard Lyons
Com esses coeficientes, o erro é sempre positivo. Com um pequeno ajuste, podemos torná-lo simétrico e reduzir pela metade o erro máximo.
Yves Daoust