Cスタンダードライブラリにおけるsin関数の実装

Cスタンダードライブラリにおける関数は様々な驚くべき実装が存在します。
rand()関数はどのようにして乱数を作っているのでしょうか。sin()関数はどのようにしてサインを計算しているのでしょうか。

ここではsin()の実装について考えてみましょう。

・sin関数の計算方法

そもそもsinを求める方法とは何があるでしょうか。


まず思いつくのはテイラー展開でしょうか。

テイラー展開をすることでsinは以下の式で表せます。

無限級数なのでどこかで丸めなければなりません。double型の桁数に合うように次数を決められれば良いのですが、その場合次数はxの値に依存することになります。
つまりxが大きいほど収束が遅くなるため、展開次数を大きくとる必要があります。



もう一つ、sinを求める方法は三角比表でしょう。

例えば以上の三角比表が与えられた場合、
sin2°はsin0°とsin5°の線形補間で計算できます。

しかし、上の表を見れば分かるように、5°刻みだと非常に粗い補間になってしまいます。
double型に使うにはもっと大きな表が必要になるでしょう。ただしそうすると、より多くのメモリを要求することになります。

また、線形補間ではなく二次補間をすることでより精度をあげることが出来るでしょう。しかしこの場合は計算時間が大きくなります。



さて、いずれの方法にしても精度と計算資源(メモリ・時間)のトレードオフの関係があると言えます。

sinの実装

では、実際の実装はどのようにしてsinを計算しているのでしょうか。

ここではIBM Libraryの実装を見てみましょう。
この実装はsinの真の値に最も近いdoubleの値を返します。

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001-2014 Free Software Foundation, Inc.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU  Lesser General Public License
 * along with this program; if not, see <http://www.gnu.org/licenses/>.
 */

//..中略..//

double
SECTION
__sin (double x)
{
  double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs, xn, a, da, db, eps, xn1,
    xn2;
  mynumber u, v;
  int4 k, m, n;
  double retval = 0;
  SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
  u.x = x;
  m = u.i[HIGH_HALF];
  k = 0x7fffffff & m;           /* no sign           */
  if (k < 0x3e500000)           /* if x->0 =>sin(x)=x */
    retval = x;
 /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
  else if (k < 0x3fd00000)
    {
      xx = x * x;
      /* Taylor series.  */
      t = POLYNOMIAL (xx) * (xx * x);
      res = x + t;
      cor = (x - res) + t;
      retval = (res == res + 1.07 * cor) ? res : slow (x);
    }                           /*  else  if (k < 0x3fd00000)    */
/*---------------------------- 0.25<|x|< 0.855469---------------------- */
  else if (k < 0x3feb6000)
    {
      u.x = (m > 0) ? big + x : big - x;
      y = (m > 0) ? x - (u.x - big) : x + (u.x - big);
      xx = y * y;
      s = y + y * xx * (sn3 + xx * sn5);
      c = xx * (cs2 + xx * (cs4 + xx * cs6));
      SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
      if (m <= 0)
        {
          sn = -sn;
          ssn = -ssn;
        }
      cor = (ssn + s * ccs - sn * c) + cs * s;
      res = sn + cor;
      cor = (sn - res) + cor;
      retval = (res == res + 1.096 * cor) ? res : slow1 (x);
    }                           /*   else  if (k < 0x3feb6000)    */
/*----------------------- 0.855469  <|x|<2.426265  ----------------------*/
  else if (k < 0x400368fd)
    {

 //..後略..//

sourceware.orgで全コードを見ることが出来ます。

まず面白い最適化として、xの大きさに応じて場合分けをしていることがあげられます。

xが小さい場合はテイラー展開を、45°程度の場合は三角比表を、などなどそれぞれの場合に異なる処理をしていることが分かります。

xが小さい場合はテイラー展開が非常に早く収束するのでテイラー展開を採用するなど、非常に高い水準でのトレードオフが見られます。


Cという言語はこのような極限の最適化を重ねて書かれた言語です。
坐臥して最速という訳ではないということでしょう。