-1

arctan（x）にMaclaurin級数を使用していますが、正しい答えが得られません。ラジアンで計算しています。これまでの関数は次のとおりです。

fp32 t32rArcTangent(fp32 number)
{
fp32 a, b, c, d;    /* Temp Variables */
fp32 t;             /* Number Temp */
uint32 i;           /* Loop Counter */

/* Time Savers */
if (b32fpcomp(number, MM_FP8INFINITY)) return((fp32)MM_PI / 2);
if (b32fpcomp(number, -MM_FP8INFINITY)) return(-(fp32)MM_PI / 2);

/* Setup */
a = 0;
b = 0;
c = 1;
d = number;
t = number * number;

/* Calculation Loop */
for (i = 0; i < MMPRVT_FP32_TRIG_LIMIT; i++)
  {
    b += d;
    if (b32fpcomp(a, b)) break;
    a = b;
    c += 2;
    d *= -1 * t / c;
  }
#ifdef DEBUG
printf("Loops: %lu\n", i);
#endif

/* Result */
return(a);

fp32 = typedef'd float

uint32 = typedef'd unsigned long int

MM_FP8INFINITYは、fp32データ型に含めることができる最大の数値です。

MM_PIは、約50桁のPI出力です。

MMPRVT_FP32_TRIG_LIMITは、結果の計算に使用できるループの最大数です。これは、何らかの理由で級数が収束しない場合に、級数の展開が無限ループになるのを防ぐためです。

これらは私が得ている結果です：

Testing arctangent(x) function.
Loops: 0
arctan(0):      0
Loops: 8
arctan(1):      0.724778414
Loops: 13
arctan(R3):     0.709577262
Loops: 6
arctan(1/R3):   0.517280579

R3は、1.732050808である3の平方根です。

これで、アークタン級数の収束半径が|x|であることがわかりました。<= 1なので、どういうわけか入力を減らす必要があると思います。問題は、arctanの場合、関数の定義域が（-INF、+ INF）であるということです。では、どうやってそれを減らすのですか？これはラジアン角度に計算されています。

それを指摘してくれてありがとう。この問題は修正され、入力の削減も行われました。これが完成して修正された関数で、これで正しい答えが得られます。

fp32 t32rArcTangent(fp32 number)
{
fp32 a, b, c, d;    /* Temp Variables */
fp32 t;             /* Number Temp */
uint32 i;           /* Loop Counter */
uint8 fr;           /* Reduction Flag */

/* Time Savers */
if (b32isInf(number) == -1) return(-(fp32)MM_PI / 2);
if (b32isInf(number) == 1) return((fp32)MM_PI / 2);
if (b32isNaN(number)) return(number);
if (b32fpcomp(number, MM_FP8INFINITY)) return((fp32)MM_PI / 2);
if (b32fpcomp(number, -MM_FP8INFINITY)) return(-(fp32)MM_PI / 2);
if (b32fpcomp(number, ONE)) return((fp32)MM_PI / 4);
if (b32fpcomp(number, -ONE)) return(-(fp32)MM_PI / 4);

/* Reduce Input */
if (number > ONE)
    {
      number = 1 / number;
      fr = 1;
    }
  else fr = 0;

/* Setup */
a = 0;
b = 0;
c = 1;
d = number;
t = number * number;

/* Calculation Loop */
for (i = 0; i < MMPRVT_FP32_TRIG_LIMIT; i++)
  {
    b += d / c;
    if (b32fpcomp(a, b)) break;
    a = b;
    c += 2;
    d *= -1 * t;
    #ifdef DEBUG
    printf("a=%g b=%g, c=%g d=%g\n", a, b, c, d);
    #endif
  }
#ifdef DEBUG
printf("Loops: %lu\n", i);
#endif

/* Result */
if (fr != 0) a = ((fp32)MM_PI / 2) - a;
return(a);
}
4

1 に答える 1

5

による除算の結果として、各ループの項に何が起こるか考えてみてくださいc。

c += 2;
d *= -1 * t / c;

最初に 1 で割っています [暗黙のうちに、この前に]、次に 3 で、次に 5 で割っています。除数。IOW、代わりに

x - 1/3*x^3 + 1/5*x^5 - 1/7*x^7 + 1/9*x^9

あなたが望むもの、あなたは計算しています

x - 1/(1*3)*x^3 + 1/(1*3*5)*x^5 - 1/(1*3*5*7)*x^7 + 1/(1*3*5*7*9)*x^9

トリックは引き続き使用できますd *= -tが、分割を移動する必要があります。

于 2012-06-17T01:46:06.903 に答える