JP2006059084A - データプロセッサ、データ処理方法及び演算処理プログラム - Google Patents

Abstract

【課題】 逆余弦関数ａｃｏｓ（ｘ）に対しては１近傍の値を高精度に得る。
【解決手段】 データプロセッサ（１）は、データ処理ユニット（２０）を有し、演算制御プログラムを実行するデータ処理ユニットは、前記演算制御プログラムに従って、入力ｘに対して関数ａｃｏｓ（ｘ）の値の計算を、多項式Ｆ（ｘ）と定数Ｃを用いてａｃｏｓ（ｘ）＝Ｆ（１−ｘ）＊ｓｑｒｔ（Ｃ＊（１−ｘ））で行うとき、Ｆ（１−ｘ）の多項式計算をホーナー法で行い、その最後の加算を留保して前記多項式Ｆ（ｘ）の定数項の値ｃｓｔと加算項の値αとを分けてストアし、更に、ｓｑｒｔ（Ｃ＊（１−ｘ））をニュートン法を用いて被加算項の値ｓｑ１と加算項の値εとに分けてストアし、（ｃｓｔ+α）＊（ｓｑ１＋ε）の展開式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εの全部又は先頭から一部を用いて、ｘの入力に対する関数ａｃｏｓ（ｘ）の値を求める演算を行う。
【選択図】 図１

Description

本発明は、逆余弦関数ａｓｏｃ（ｘ）の値を求めるための演算制御技術に関し、例えばマイクロコンピュータに適用して有効な技術に関する。

ｓｉｎ（ｘ）、ｃｏｓ（ｘ）、ｅｘｐ（ｘ）、ｌｏｇ（ｘ）などの数学関数は応用分野によっては多用される。それら関数の値を演算するプログラムはＣなど高級言語の標準ライブラリに含まれている。そのようなプログラムでは高速かつ高精度な近似関数による演算制御を行うことができる。近似関数を作る場合、近似区間を定め、その区間内で目標関数を精度よく近似する多項式や有理式を決め、その近似式で値を計算する。また近似区間外はその近似区間内へ還元して値を計算する。上記ｓｉｎ（ｘ）などに対しては容易に多項式近似関数などを導出でき、また、その近似関数で得られる関数値は真の値に極めて近いものとなる。しかしながら、逆余弦関数ａｃｏｓ（ｘ）に対しては１近傍の近似式を作るのが難しい。その関数は、１近傍で正の値が急激に０に落ちていく関数形状であり、多項式で近似するのにはなじみ難い。多項近似について記載された文献として例えば非特許文献１がある。

CODY, W. J., AND WAITE, W. Software Manual for the Elementary Functions. Prentice-Hall, Englewood Cliffs, N.J., 1980.

本発明者は逆余弦関数ａｃｏｓ（ｘ）に対しては１近傍の近似式について検討した。１近傍におけるａｃｏｓの関数形状は平方根の関数ｓｑｒｔ（ｘ）に近似している。Ｆ（ｙ）＝ａｃｏｓ（１−ｙ）／ｓｑｒｔ（ｙ）とおく。Ｆ（ｙ）は多項近似することができる。ｓｑｒｔ（ｙ）の関数演算命令は多くのプロセッサの命令セットに含まれているから、その値は容易に計算できる。そうすると、ａｃｏｓ（１−ｙ）＝Ｆ（ｙ）＊ｓｑｒｔ（ｙ）の演算によってａｃｏｓ（１−ｙ）の値を演算できることになる。＊は乗算記号である。

しかしながら、単にそれだけでは演算精度に問題のあることが本発明者によって明らかにされた。すなわち、計算を出力のデータ型以上の精度で行い、最後に出力のデータ型に丸めるというやり方を行えば精度上問題はないのであるが、出力のデータ型で計算しなければならない（または拡張精度の演算命令はない）という制約があると、Ｆ（ｙ）の計算とｓｑｒｔ（ｙ）の計算で各々丸め誤差が発生し、その後の乗算で更に丸め誤差が発生するからである。

本発明の目的は、逆余弦関数ａｃｏｓ（ｘ）に対しては１近傍の値を高精度に得ることができる演算制御技術を提供することにある。

本発明の前記並びにその他の目的と新規な特徴は本明細書の記述及び添付図面から明らかになるであろう。

本願において開示される発明のうち代表的なものの概要を簡単に説明すれば下記の通りである。

〔１〕データプロセッサは、データ処理ユニットと、前記データ処理ユニットが実行する演算制御プログラムを保有するプログラムメモリと、前記演算制御プログラムにしたがって前記データ処理ユニットがアクセスするワークメモリと、を有する。前記データ処理ユニットは、前記演算制御プログラムに従って、入力ｘに対して関数ａｃｏｓ（ｘ）の値の計算を、多項式Ｆ（ｘ）と定数Ｃを用いてａｃｏｓ（ｘ）＝Ｆ（１−ｘ）＊ｓｑｒｔ（Ｃ＊（１−ｘ））で行うとき、Ｆ（１−ｘ）の多項式計算をホーナー法で行い、その最後の加算を留保して前記多項式Ｆ（ｘ）の定数項の値ｃｓｔと加算項の値αとを分けてストアし、更に、ｓｑｒｔ（Ｃ＊（１−ｘ））をニュートン法を用いて被加算項の値ｓｑ１と加算項の値εとに分けてストアし、（ｃｓｔ+α）＊（ｓｑ１＋ε）の展開式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εの全部又は先頭から一部を用いて、ｘの入力に対する関数ａｃｏｓ（ｘ）の値を求める演算を行う。要するに、Ｆ（１−ｙ）＝ｃｓｔ+α、ｓｑｒｔ（Ｃ＊（１−ｘ））＝ｓｑ１＋εのように、Ｆ（１−ｙ）とｓｑｒｔ（Ｃ＊（１−ｘ））を各々２数の和とする。ａｃｏｓ（ｘ）の１近傍における関数形状は特にＣ＝２の平方根の関数ｓｑｒｔ（Ｃ＊（１−ｘ））に近似しているから、Ｃ＝２とすれば、ｃｓｔ＝１．０として差し支えない。αはホーナー法により、εはニュートン法により、夫々計算精度の高い比較的小さな値とされる。このとき、上記恒等式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εについて見ると、式の４項は左側から大きい順に並んでいる。この式を計算してａｃｏｓ（ｘ）の値を求める。このとき、右側の方から加算していく。そうすると、丸めの効果が次の加算に反映されず、最終結果は高精度が期待される。尚、式の最後の項α＊εは値が極端に小さく、加算しなくてもなんら支障はない。

本発明の具体的な形態として、ホーナー法によるＦ（１−ｘ）の多項式計算において、多項式をＦ（ｙ）＝ｃ０＋ｃ１＊ｙ＋ｃ２＊ｙ＾２＋ｃ３＊ｙ＾３…とすると、前記多項式を、Ｆ（ｙ）＝ｃ０＋ｙ＊（ｃ１＋ｙ＊（ｃ２＋ｙ＊（ｃ３…）））とし、入力ｙに対する前記多項式の値をｙの高次側よりホーナー法により演算し、保留する最後の積和演算をｃ０＋ｙ＊ａｎｓとし、ｃ０＝ｃｓｔ、ｙ＊ａｎｓ=αとする。ｃ０には演算誤差はない。

また、前記データ処理ユニットは平方根の関数ｓｑｒｔの演算処理にてｓｑｒｔ（Ｃ＊ｙ）の値ｓｑを演算し、その値ｓｑに対してニュートン法を適用し、ｓｑ＝０．５＊（ｓｑ＋Ｃ＊ｙ／ｓｑ）をｓｑ＝ｓｑ＋０．５＊（（Ｃ＊ｙ−ｓｑ＊ｓｑ）ｓｑ）と変形し、ｓｑ１＝ｓｑ、ε＝０．５＊（（Ｃ＊ｙ−ｓｑ＊ｓｑ）ｓｑ）とする。

本発明の更に具体的な形態として、前記定数Ｃが２であり、前記多項式Ｆ（ｘ）の値ｃｓｔは１．０である。前記定数Ｃを１とするときは、前記多項式Ｆ（ｘ）の定数係数を相対的に大きな値ｒｏｏｔ２Ｋと相対的に小さい値ｒｏｏｔ２ｆに分け、前記値ｃｓｔとして前記値ｒｏｏｔ２Ｋを採用し、前記値αに前記値ｒｏｏｔ２ｆを含める。

〔２〕データ処理方法は、コンピュータ装置が演算制御プログラムを実行することにより、入力ｘに対して関数ａｃｏｓ（ｘ）の値を、多項式Ｆ（ｘ）と定数Ｃを用いてａｃｏｓ（ｘ）＝Ｆ（１−ｘ）＊ｓｑｒｔ（Ｃ＊（１−ｘ））の計算により取得する方法であって、Ｆ（１−ｘ）の多項式計算をホーナー法で行い、その最後の加算を留保して前記多項式Ｆ（ｘ）の定数項の値ｃｓｔと加算項の値αとを分けてストアする処理と、ｓｑｒｔ（Ｃ＊（１−ｘ））をニュートン法を用いて被加算項の値ｓｑ１と加算項の値εとに分けてストアする処理と、（ｃｓｔ+α）＊（ｓｑ１＋ε）の展開式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εの全部又は先頭から一部を用いて、ｘの入力に対する関数ａｃｏｓ（ｘ）の値を演算する処理と、を含む。

〔３〕コンピュータ装置によって実行される演算制御プログラムは、入力ｘに対して関数ａｃｏｓ（ｘ）の値のを、多項式Ｆ（ｘ）と定数Ｃを用いてａｃｏｓ（ｘ）＝Ｆ（１−ｘ）＊ｓｑｒｔ（Ｃ＊（１−ｘ））の計算により取得するデータ処理の制御記述を有し、前記データ処理は、Ｆ（１−ｘ）の多項式計算をホーナー法で行い、その最後の加算を留保して前記多項式Ｆ（ｘ）の定数項の値ｃｓｔと加算項の値αとを分けてストアする処理と、ｓｑｒｔ（Ｃ＊（１−ｘ））をニュートン法を用いて被加算項の値ｓｑ１と加算項の値εとに分けてストアする処理と、（ｃｓｔ+α）＊（ｓｑ１＋ε）の展開式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εの全部又は先頭から一部を用いて、ｘの入力に対する関数ａｃｏｓ（ｘ）の値を演算する処理と、を含む。

本願において開示される発明のうち代表的なものによって得られる効果を簡単に説明すれば下記の通りである。

すなわち、逆余弦関数ａｃｏｓ（ｘ）に対しては１近傍の値を高精度に得ることができる。

図１には本発明の一例に係るマイクロコンピュータが示される。同図に示されるマイクロコンピュータ１は、特に制限されないが、公知の半導体集積回路製造技術によって単結晶シリコンのような１個の半導体基板（半導体チップ）に形成される。このマイクロコンピュータ１は、データ処理ユニット２０として、例えば、整数演算を行なう中央処理装置（ＣＰＵ）２と共に、浮動小数点演算を行なう浮動小数点演算装置（ＦＰＵ）３を有する。更に、ディジタル信号処理で多用される積和演算に特化した積和演算装置４を有する。

中央処理装置２は内部アドレスバス１４及び内部データバス１３に結合される。中央処理装置２は、特に制限されないが、汎用レジスタや算術論理演算器で代表される演算部２Ａと、プログラムカウンタなどの制御用レジスタ群、そして命令のフェッチや解読並びに命令実行手順を制御したり演算制御を行う命令制御部２Ｂなどを有する。前記内部バス１３，１４には、前記ＣＰＵ２が実行する演算制御プログラムなどを保有するプログラムメモリ（ＰＧＭＭ）２１と、前記演算制御プログラムにしたがって前記ＣＰＵ２がアクセスするワークメモリ（ＷＲＫＭ）２２とを有する。中央処理装置２は、内部バス１３，１４に接続されたプログラムメモリ２１から命令をフェッチし、その命令を解読し、解読結果に応ずる制御信号を生成することにより、当該命令に応じたデータ処理を行う。

浮動小数点演算装置（ＦＰＵ）３及び積和演算装置４は内部データバス１３に結合される。浮動小数点演算装置（ＦＰＵ）３及び積和演算装置４は図示を省略する演算回路と共にデータレジスタを有し、このデータレジスタにメモリから演算データがロードされ、演算結果データは、そのレジスタからメモリにストアされる。前記ロード、ストアなどのためのアドレッシング動作はＣＰＵ２が行なう。したがって、ＦＰＵ３及び積和演算装置４はメモリアクセスのためのメモリアドレシング能力を備える必要はない。これは、ＦＰＵ３及び積和演算装置４によるメモリアドレシング回路の必要性を取り除いてチップ面積を節約するためである。

内部データバス１３及び内部アドレスバス１４はバスステートコントローラ６に結合される。マイクロコンピュータ１による外部アクセスは、前記バスステートコントローラ６に接続された外部バスインタフェース回路８で行う。外部バスインタフェース回路８は外部データバス１８及び外部アドレスバス１７に接続される。また、前記バスステートコントローラ６には、周辺データバス１６及び周辺アドレスバス１５を介して、例えば、クロック発生回路７、システムコントローラ１２、シリアルコミュニケーションインタフェースコントローラ（ＳＣＩ）９、タイマ１０及びＡ／Ｄコンバータ１１が結合される。それら周辺回路はデータレジスタや制御レジスタを有し、それらレジスタは前記バスステートコントローラ６を介してＣＰＵ２によってアクセスされる。

ＣＰＵ２が管理する内部メモリ空間、外部メモリ空間、前記ＳＣＩ９，タイマ１０，Ａ／Ｄコンバータ１１などの周辺回路に対するアドレスエリアの割り当ては予め決定されている。前記バスステートコントローラ６は、アクセスエリア毎にアクセスサイクル数やバス幅などがＣＰＵ２によって設定される図示を省略するバスコントロールレジスタを有し、ＣＰＵ２からのアクセスアドレスで指定されるメモリエリアに対するバス幅やアクセスサイクル数などのバス制御を行なって、バスサイクルを起動する。ＣＰＵ２からのアクセスの指示はバスコマンド２３としてバスステートコントローラ６に与えられる。バスコマンド２３にはアクセスサイズの指定や、リード、ライト、メモリアクセス等のストローブ信号が含まれる。

割込みコントローラ５は複数の割込み要求に対する優先制御やマスク制御などを行なって割込み信号２５をＣＰＵ２に与える。割込み要求は、バスステートコントローラ６からのバスエラー信号２４のほか、ＳＣＩ９、タイマ１０、Ａ／Ｄコンバータ１１そして外部からの図示を省略する割込み要求信号によって与えられる。

前記マイクロコンピュータ１は、特に制限されないが、クロック発生回路７から出力されるクロック信号ＣＬＫ０〜ＣＬＫ２に同期動作される。前記システムコントローラ１２はマイクロコンピュータ１の内蔵モジュールに対する動作停止の制御などを行う。

前記データ処理ユニット２０は、プログラムメモリ２１が保有する前記演算制御プログラムに従って、入力ｘに対して関数ａｃｏｓ（ｘ）の値の計算を、多項式Ｆ（ｘ）と定数Ｃを用いてａｃｏｓ（ｘ）＝Ｆ（１−ｘ）＊ｓｑｒｔ（Ｃ＊（１−ｘ））で行うとき、Ｆ（１−ｘ）の多項式計算をホーナー法で行い、その最後の加算を留保して前記多項式Ｆ（ｘ）の定数項の値ｃｓｔと加算項の値αとを分けてストアし、更に、ｓｑｒｔ（Ｃ＊（１−ｘ））をニュートン法を用いて被加算項の値ｓｑ１と加算項の値εとに分けてストアし、（ｃｓｔ+α）＊（ｓｑ１＋ε）の展開式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εの全部又は先頭から一部を用いて、ｘの入力に対する関数ａｃｏｓ（ｘ）の値を求める演算を行う。要するに、Ｆ（１−ｘ）＝ｃｓｔ+α、ｓｑｒｔ（Ｃ＊（１−ｘ））＝ｓｑ１＋εのように、Ｆ（１−ｘ）とｓｑｒｔ（Ｃ＊（１−ｘ））を各々２数の和とする。ａｃｏｓ（ｘ）の１近傍における関数形状は特にＣ＝２の平方根の関数ｓｑｒｔ（Ｃ＊（１−ｘ））に近似しているから、Ｃ＝２とすれば、ｃｓｔ＝１．０として差し支えない。αはホーナー法により、εはニュートン法により、夫々計算精度の高い比較的小さな値とされる。このとき、上記恒等式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εについて見ると、式の４項は左側から大きい順に並んでいる。この式を計算してａｃｏｓ（ｘ）の値を求める。このとき、右側の方から加算していく。そうすると、丸めの効果が次の加算に反映されず、最終結果は高精度が期待される。尚、式の最後の項α＊εは値が極端に小さく、加算しなくてもなんら支障はない。

以下に上記演算制御方法を説明する。ここでは、演算速度を確保しつつ、演算精度を大幅に向上する方法を提案する。単精度浮動小数点数で値を返すａｃｏｓ（ｘ）を単精度浮動小数点演算で計算する例を説明するが、倍精度浮動小数点数で値を返すａｃｏｓ（ｘ）を倍精度浮動小数点演算で計算する場合にもそのまま適用することができる。

まず、基本の発想を示すために
Ａ＝（１＋α）＊（ｓｑ＋ε）…式１
の式を考える。α、ｓｑ、εは正の値である。式1の右辺を展開すると、
Ａ＝ｓｑ＋α＊ｓｑ＋ε＋α＊ε…式２
のようになる。ここで、αとｓｑは１よりも１桁以上小さく、εはｓｑよりも８桁以上小さいとする。値は総て単精度浮動小数点数に格納されるものとする。値Ａを式１の通りに計算すると、２つの加算で丸め誤差が発生し、最後の乗算でさらに丸め誤差が発生する。しかし、式２で計算し、加算を右側から行うとすれば、誤差の発生は最小限に抑えられる。式２の４項は大きい順に配置されていて、左側２つの項の比がかなり大きいからである。さらに最右翼の加算は省略できるであろう。最右翼の項は、最終のＡの値に寄与できない程度に小さいからである。

以上を基本の発想として、ａｃｏｓ（ｘ）の高精度算法を考える。近似すべき目標関数を、
Ｆ（ｙ）＝ａｃｏｓ（１−ｙ）／ｓｑｒｔ（Ｃ＊ｙ）…式３
とする。Ｃ＝２である。本明細書においてｙ＝１−ｘとする。式３の分子のｓｑｒｔ（）の中のＣを２とすると、分母と分子の関数形状は殆ど等しくなる。要するに、ａｃｏｓ（１−ｙ）に原点位置で接する放物線がｓｑｒ（２＊ｙ）となる。従って、Ｆ（ｙ）の定数係数を１．０に固定した近似多項式を導出できる。例えばＦ（ｙ）＝ｃ０＋ｃ１＊ｙ＋ｃ２＊ｙ＾２＋ｃ３＊ｙ＾３…とすると、前記近似多項式は図２の係数を持つ。近似区間は［０〜０．３］である。図２において最初の値は限界誤差（絶対）、残りが係数で逆順に表示されている。最下段が定数係数である。

図３はＦ（ｙ）を倍精度で計算したときの誤差グラフを示し、図４はａｃｏｓ（１−ｘ）を倍精度で計算したときの誤差グラフを示す。表示区間は近似区間と同じで、［０〜０．３］である。共に絶対誤差で表示している。

さて、Ｆ（ｙ）をホーナー法で計算する場合、Ｆ（ｙ）＝ｃ０＋ｙ＊（ｃ１＋ｙ＊（ｃ２＋ｙ＊（ｃ３＋ｙ＊（ｃ４＋ｙ＊（ｃ５）））））とすると、このＦ（ｙ）を素直にホーナー法で計算すれば以下に示す演算、
ａｎｓ＝ｃ５、
ａｎｓ＝ｃ４＋ａｎｓ＊ｙ、
ａｎｓ＝ｃ３＋ａｎｓ＊ｙ、
ａｎｓ＝ｃ２＋ａｎｓ＊ｙ、
ａｎｓ＝ｃ１＋ａｎｓ＊ｙ、
ａｎｓ＝１．０＋ａｎｓ＊ｙ、
が行われる。この最後の積和の計算（ａｎｓ＝１．０＋ａｎｓ＊ｙ）を留保する。代わりにα=ａｎｓ＊ｙを求めておく。このαは上記式１、式２のαに丁度対応している。これにより、式２で計算するａｃｏｓ（ｘ）の半分の準備が整う。

次に、ｓｑ＝ｓｑｒｔ（２＊ｙ）の計算を行う。このｓｑｒｔ（２＊ｙ）はＩＥＥＥＥの規格通りの計算をするものとする。結果の値は単精度浮動小数点の値として最善のものになる。しかし、それ以上にｓｑの精度を確保したい。そこでニュートン法を適用して精度を上げる。すなわち、
ｓｑ＝０．５＊（ｓｑ＋２＊ｙ／ｓｑ）
とし、これをｓｑ＝ｓｑ＋０．５＊（（２＊ｙ−ｓｑ＊ｓｑ）ｓｑ）と変形する。変形前の算法であれば新しいｓｑは前のｓｑと同じ値となり無意味である。なぜ同じ値かといえば、以前のｓｑが最善の近似値だったからである。一方、変形後の算法で第２項目の加算を留保すれば精度がほぼ２倍になる。つまり、ε＝（２＊ｙ―ｓｑ＊ｓｑ）／（２＊ｓｑ）とｓｑの２つで新しいｓｑを考えるということである。なお、上記変形式において、（２＊ｙ−ｓｑ＊ｓｑ）は２＊ｙを演算した後で（２＊ｙ）−ｓｑ＊ｓｑを積和命令で計算するものとし、積和命令はｓｑ＊ｓｑの積の結果を総て和（実際は差）の演算に参加させるという前提が必要である。

これで、式２で計算するａｃｏｓ（ｘ）の残りの半分の準備が整ったことになり、これらを用いて式２の演算を行う。

図５にはａｃｏｓ（ｘ）の値を求めるための上記演算を採用したプログラム記述の例が示される。図５の記法はＣ言語に準拠し、ｆｌｏａｔは単精度浮動小数点数を示す。図６には図５と同一処理をＣ言語による実際の記述法に則した別の記法で示す。図５、６においてｙ２＝ｙ＋ｙは、前記式３のｓｑｒｔ（Ｃ＊ｙ）においてＣ＝２としたとき、２＊ｙを便宜上ｙ＋ｙ（＝ｙ２）として表現したものである。図７にはホーナー法における最後の加算演算の留保とニュートン法による演算を行わない場合の比較例としてのプログラム記述の例が示される。図８には図７と同一処理をＣ言語による実際の記述法に則した別の記法で示す。

図５、６に代表される演算方法による演算結果に対する演算精度の評価結果について説明する。図５、６に代表される演算方法をｎｅｗａｃｏｓ法、図７、８に代表される比較例にかかる演算方法をｏｌｄａｃｏｓ法と称する。評価を行うのに、近似区間［０．７〜１．０］において０．００２刻みにサンプル点ｘを定め、ｎｅｗ＿ａｃｏｓ（ｘ）−（ｆｌｏａｔ）ａｃｏｓ（ｘ）のような差分を調べた。すなわち、倍精度の浮動小数点数の値を返すライブラリ関数ａｃｏｓ（ｘ）が返した値を単精度浮動小数点数に丸めたしたもの（（ｆｌｏａｔ）ａｃｏｓ（ｘ））を基準値とし、ｎｅｗａｃｏｓ法による関数値（ｎｅｗ＿ａｃｏｓ（ｘ））と差分を取った。その結果は図９に示される。ｏｌｄａｃｏｓ法による関数値（ｏｌｄ＿ａｃｏｓ（ｘ））と差分も取った。その結果は図１０に示される。ｏｌｄａｃｏｓ法は、１／３のサンプル点で値が１つずれている。一方ｎｅｗａｃｏｓ法は総てのサンプル点で基準の（ｆｌｏａｔ）ａｃｏｓ（ｘ）と同じであった。［０．８〜１．０］を０．００１刻みとした別のサンプル点群で調べた結果を図１１と図１２に示すが、この場合も上記と別のサンプル点群で調べても、同様な傾向であった。図１１、図１２の結果は図９、図１０の結果に比べて、誤差が少し増えているｏｌｄａｃｏｓ（ｘ）では、値が２つずれているサンプル点が現れている。また、ｎｅｗａｃｏｓ（ｘ）では、値が１つずれているものが現れている。しかし基本的には図９、図１０と同様の傾向といえる。いずれにしても、式２に従って精度確保を目指したｎｅｗａｃｏｓ法による演算制御によれば極めて良好な演算精度を得ることができる。誤差の混入を最小限におさえた演算制御方法だからである。

次に演算速度について考察する。前述の図７及び図８のｏｌｄ＿ａｃｏｓｆの関数は前記式１にしたがって計算するコードで、単にホーナー法で計算しているだけである。図５及び図６のｎｅｗ＿ａｃｏｓｆの関数が前記式２に従って演算を行うためのコードである。２つの関数、値の定義／使用の主系列部分を普通の位置に、そうでない系列部分は右にずらせて配置している。入力はｙである。ｓｑは平方根演算で求まるが、平方根演算命令はあり、その完了クロック数は１２とする。積和命令や乗算命令の完了クロック数は３とする。関数ｏｌｄ＿ａｃｏｓではホーナー法の計算部分が演算主系列を構成する。前の積和演算の結果を次の積和演算で使用するリカレンスの演算構造である。上の関数の演算部分のクロック数は１８（＝６＊３）となる。

関数ｎｅｗ＿ａｃｏｓでは、平方根演算とその後の除算部分が演算主系列に含まれる。そして除算の完了を１２クロックとすれば、下の関数の演算部分のクロック数は３３（＝２＊１２＋３＊３）となる。

従って、演算主系列の意味でのクロック数は１５違うことになる。実際に２つのコードをある２並列のスーパースカラプロセッサで実行させ、関数としてのクロック数を計ると、２５クロックと４０クロックであった。つまり机上計算の差と丁度同じになった。

以上の説明では、ａｃｏｓ＝Ｆ（１−ｘ）／ｓｑｒｔ（２＊（１−ｘ））として説明した。この場合多項式Ｆ（ｘ）の定数係数を１．０とできた。しかしｓｑｒｔにおけるＣは２である必要はない。例えば１でもいい。この場合多項式Ｆ（ｘ）の定数係数はｓｑｒｔ（２）に近い値となる。これを相対的に大きい値ｒｏｏｔ２ｋと相対的に小さいｒｏｏｔ２の２つの和で管理するとし、ホーナー法の計算を
ａｎｓ＝ｃ５
ａｎｓ＝ｃ４＋ａｎｓ＊ｙ
ａｎｓ＝ｃ３＋ａｎｓ＊ｙ
ａｎｓ＝ｃ２＋ａｎｓ＊ｙ
ａｎｓ＝ｃ１＋ａｎｓ＊ｙ
α＝ｒｏｏｔ２ｆ＋ａｎｓ＊ｙ
ａｎｓ＝ｒｏｏｔ２ｋ+α
とし、最後の演算を留保する。更に式１の関係を
Ａ＝（ｒｏｏｔ２ｋ＋α）＊（ｓｑ+ε）
のように変更し、対応する式２の関係を以下の
Ａ＝ｒｏｏｔ２ｋ＊ｓｑ＋α＊ｓｑ＋ｒｏｏｔ２ｋ＊ε＋α＊ε
のように変形する。そして、この式を右側から加算していく。勿論、最右の項は無視してもいい。

以上本発明者によってなされた発明を実施形態に基づいて具体的に説明したが、本発明はそれに限定されるものではなく、その要旨を逸脱しない範囲において種々変更可能であることは言うまでもない。

例えば、データプロセッサは１チップのマイクロコンピュータに限定されず、マルチチップのプロセッサであってもよい。したがって、個別チップ化されたプログラムメモリ、ワークメモリ及び演算処理ユニットを回路基板に搭載してデータプロセッサを構成してもよい。

本発明の一例に係るマイクロコンピュータのブロック図である。 Ｆ（ｙ）＝ｃ０＋ｃ１＊ｙ＋ｃ２＊ｙ＾２＋ｃ３＊ｙ＾３…の係数を例示する説明図である。 Ｆ（ｙ）を倍精度で計算したときの誤差を示す特性図である。 ａｃｏｓ（１−ｘ）を倍精度で計算したときの誤差を示す特性図である。 ｎｅｗａｃｏｓ法を採用した演算制御プログラムの記述例を示す説明図である。 図５と同一処理をＣ言語による実際の記述法に則して示した説明図である。 比較例としてのｏｌｄａｃｏｓ法による演算制御プログラムの記述例を示す説明図である。 図７と同一処理をＣ言語による実際の記述法に則して示した説明図である。 倍精度の浮動小数点数の値を返すライブラリ関数ａｃｏｓ（ｘ）が返した値を単精度浮動小数点数に丸めた基準値に対するｎｅｗａｃｏｓ法による関数値との差分を誤差として示す誤差特性図である。 図９と同様の基準値に対するｏｌｄａｃｏｓ法による関数値と差分を誤差として示す誤差特性図である。 別のサンプル点群に対し図９と同様の基準値に対してｎｅｗａｃｏｓ法による関数値との差分を誤差として示す誤差特性図である。 別のサンプル点群に対し図９と同様の基準値に対してｏｌｄａｃｏｓ法による関数値との差分を誤差として示す誤差特性図である。

符号の説明

１ マイクロコンピュータ
２ 中央処理装置
３ 浮動小数点演算装置
４ 積和演算装置
２０ データ処理ユニット
２１ プログラムメモリ
２２ ワークメモリ

Claims (11)

1. データ処理ユニットと、前記データ処理ユニットが実行する演算制御プログラムを保有するプログラムメモリと、前記演算制御プログラムにしたがって前記データ処理ユニットがアクセスするワークメモリと、を有し、
   前記データ処理ユニットは、前記演算制御プログラムに従って、入力ｘに対して関数ａｃｏｓ（ｘ）の値の計算を、多項式Ｆ（ｘ）と定数Ｃを用いてａｃｏｓ（ｘ）＝Ｆ（１−ｘ）＊ｓｑｒｔ（Ｃ＊（１−ｘ））で行うとき（但し＊は乗算記号）、Ｆ（１−ｘ）の多項式計算をホーナー法で行い、その最後の加算を留保して前記多項式Ｆ（ｘ）の定数項の値ｃｓｔと加算項の値αとを分けてストアし、更に、ｓｑｒｔ（Ｃ＊（１−ｘ））をニュートン法を用いて被加算項の値ｓｑ１と加算項の値εとに分けてストアし、（ｃｓｔ+α）＊（ｓｑ１＋ε）の展開式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εの全部又は先頭から一部を用いて、ｘの入力に対する関数ａｃｏｓ（ｘ）の値を演算するデータプロセッサ。
2. ホーナー法によるＦ（１−ｘ）の多項式計算において、多項式をＦ（ｙ）＝ｃ０＋ｃ１＊ｙ＋ｃ２＊ｙ＾２＋ｃ３＊ｙ＾３…とすると（但しｙ＝１−ｘ）、前記多項式を、Ｆ（ｙ）＝ｃ０＋ｙ＊（ｃ１＋ｙ＊（ｃ２＋ｙ＊（ｃ３…）））とし、入力ｙに対する前記多項式の値をｙの高次側よりホーナー法により演算し、保留する最後の積和演算をｃ０＋ｙ＊ａｎｓとし、ｃ０＝ｃｓｔ、ｙ＊ａｎｓ=αとする請求項１記載のデータプロセッサ。
3. 前記データ処理ユニットは平方根の関数ｓｑｒｔの演算処理にてｓｑｒｔ（Ｃ＊ｙ）の値ｓｑを演算し、その値ｓｑに対してニュートン法を適用し、ｓｑ＝０．５＊（ｓｑ＋Ｃ＊ｙ／ｓｑ）をｓｑ＝ｓｑ＋０．５＊（（Ｃ＊ｙ−ｓｑ＊ｓｑ）ｓｑ）と変形し、ｓｑ１＝ｓｑ、ε＝０．５＊（（Ｃ＊ｙ−ｓｑ＊ｓｑ）ｓｑ）とする請求項１記載のデータプロセッサ。
4. 前記定数Ｃが２であり、前記多項式Ｆ（ｘ）の値ｃｓｔは１．０である請求項１記載のデータプロセッサ。
5. 前記定数Ｃは１であり、前記多項式Ｆ（ｘ）の値を相対的に大きな値ｒｏｏｔ２Ｋと相対的に小さい値ｒｏｏｔ２ｆに分け、前記値ｃｓｔとして前記値ｒｏｏｔ２Ｋを採用し、前記値αに前記値ｒｏｏｔ２ｆを含める請求項１記載のデータプロセッサ。
6. コンピュータ装置が演算制御プログラムを実行することにより、入力ｘに対して関数ａｃｏｓ（ｘ）の値を、多項式Ｆ（ｘ）と定数Ｃを用いてａｃｏｓ（ｘ）＝Ｆ（１−ｘ）＊ｓｑｒｔ（Ｃ＊（１−ｘ））の計算により取得するデータ処理方法であって、
   Ｆ（１−ｘ）の多項式計算をホーナー法で行い、その最後の加算を留保して前記多項式Ｆ（ｘ）の定数項の値ｃｓｔと加算項の値αとを分けてストアする処理と、
   ｓｑｒｔ（Ｃ＊（１−ｘ））をニュートン法を用いて被加算項の値ｓｑ１と加算項の値εとに分けてストアする処理と、
   （ｃｓｔ+α）＊（ｓｑ１＋ε）の展開式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εの全部又は先頭から一部を用いて、ｘの入力に対する関数ａｃｏｓ（ｘ）の値を演算する処理と、を含むデータ処理方法。
7. 前記定数Ｃが２であり、前記多項式Ｆ（ｘ）の値ｃｓｔは１．０である請求項６記載のデータ処理方法。
8. 前記定数Ｃは１であり、前記多項式Ｆ（ｘ）の定数係数を相対的に大きな値ｒｏｏｔ２Ｋと相対的に小さい値ｒｏｏｔ２ｆに分け、前記値ｃｓｔとして前記値ｒｏｏｔ２Ｋを採用し、前記値αに前記値ｒｏｏｔ２ｆを含める請求項６記載のデータ処理方法。
9. コンピュータ装置によって実行される演算制御プログラムであって、入力ｘに対して関数ａｃｏｓ（ｘ）の値を、多項式Ｆ（ｘ）と定数Ｃを用いてａｃｏｓ（ｘ）＝Ｆ（１−ｘ）＊ｓｑｒｔ（Ｃ＊（１−ｘ））の計算により取得するデータ処理の制御記述を有し、前記データ処理は、
   Ｆ（１−ｘ）の多項式計算をホーナー法で行い、その最後の加算を留保して前記多項式Ｆ（ｘ）の定数項の値ｃｓｔと加算項の値αとを分けてストアする処理と、
   ｓｑｒｔ（Ｃ＊（１−ｘ））をニュートン法を用いて被加算項の値ｓｑ１と加算項の値εとに分けてストアする処理と、
   （ｃｓｔ+α）＊（ｓｑ１＋ε）の展開式ｃｓｔ＊ｓｑ１＋α＊ｓｑ１＋ｃｓｔ＊ε＋α＊εの全部又は先頭から一部を用いて、ｘの入力に対する関数ａｃｏｓ（ｘ）の値を演算する処理と、を含む演算処理プログラム。
10. 前記定数Ｃが２であり、前記多項式Ｆ（ｘ）の値ｃｓｔは１．０である請求項９記載の演算処理プログラム。
11. 前記定数Ｃは１であり、前記多項式Ｆ（ｘ）の定数係数を相対的に大きな値ｒｏｏｔ２Ｋと相対的に小さい値ｒｏｏｔ２ｆに分け、前記値ｃｓｔとして前記値ｒｏｏｔ２Ｋを採用し、前記値αに前記値ｒｏｏｔ２ｆを含める請求項９記載の演算処理プログラム。

Images