(57)【要約】 【課題】  演算装置による有効な演算ビット数を拡張することなく関数の近似多項式による演算精度を向上させる。【解決手段】  データ処理ユニット(２０)は、演算制御プログラムに従って所定の関数の値を演算するとき、近似多項式F(x)を分解した和の多項式F1(x)+F2(x)を用いる。前記式F1(x)は前記データ処理ユニットの演算精度に対して正確に計算可能な次数と係数を有しF1(x)＝BIG＋SMALLとして表現可能である。F1(x)の絶対値はF2(x)の絶対値よりも大きい。BIGの絶対値はSMALLの絶対値よりも大きい。前記多項式F1(x)＋F2(x)を用いた演算において、F(x)＝BIG＋SMALL＋F2(x)を求めるとき、BIG＋(SMALL＋F2(ｘ))の演算を行う。これにより、絶対値の小さいもの同士の加算で生じた丸め誤差はBIGに対して1度だけ作用する。丸め誤差が累積して大きな演算誤差を生ずる方向に最終解が丸められる可能性を低減することができる。【選択図】        図１
2006323710.tif 000002
【特許請求の範囲】
【請求項１】
データ処理ユニットと、前記データ処理ユニットが実行する演算制御プログラムを保有するプログラムメモリと、前記演算制御プログラムにしたがって前記データ処理ユニットがアクセスするワークメモリと、を有し、
前記データ処理ユニットは、前記演算制御プログラムに従って所定の関数の値を演算するとき、近似多項式F(x)を分解した和の多項式F1(x)+F2(x)を用い、
前記式F1(x)は前記データ処理ユニットの演算精度に対して正確に計算可能な次数と係数を有しF1(x)＝BIG＋SMALLとして表現可能であり、F1(x)の絶対値はF2(x)の絶対値よりも大きく、BIGの絶対値はSMALLの絶対値よりも大きく、
前記多項式F1(x)＋F2(x)を用いた演算において、F(x)＝BIG＋SMALL＋F2(x)を求めるとき、BIG＋(SMALL＋F2(ｘ))を行うデータプロセッサ。
【請求項２】
前記所定の関数の近似多項式は、F(x)＝x＋c1*x^3＋c2*x^5＋c3*x^7＋c4*x^9＋c5*x^11であり、
これに対応する和の多項式はF(x)＝(x＋C1*x^3)＋((c1−C1)*x^3＋c2*x^5＋c3*d1^7＋c4*x^9＋c5*x^11)であり、
F1(x)＝(x＋C1*x^3)であり、
F2(x)＝((c1−C1)*x^3＋c2*x^5＋c3*d1^7＋c4*x^9＋c5*x^11)であり、
C1は前記演算精度の範囲内のビット数で完全に表現できる浮動小数点数であり、
F1(x)/x＝(1＋C1*x^2)の部分演算において、1＋C1*x^2の演算結果の上位桁の値をB1に格納し、下位桁の値をS1に格納し、
BIG＝x*B1、SMALL＝x*S1とする請求項１記載のデータプロセッサ。
【請求項３】
前記所定の関数の近似多項式は、F(x)＝c0＋c1*x^2＋c2*x^4＋c3*x^6＋c4*x^8＋c5*x^10であり、
これに対応する和の多項式はF(x)＝(c0＋C1*x^2)＋((c1−C1)*x^2＋c2*x^4＋c3*d1^6＋c4*x^8＋c5*x^10)であり、
F1(x)＝(c0＋C1*x^2)であり、
F2(x)＝((c1−C1)*x^2＋c2*x^4＋c3*d1^6＋c4*x^8＋c5*x^10)であり、
C1は前記演算精度の範囲内のビット数で完全に表現できる浮動小数点数であり、
F1(x)＝(c0＋C1*x^2)の部分演算において、c0＋C1*x^2の演算結果の上位桁の値をB1に格納し、下位桁の値をS1に格納し、
BIG＝B1、SMALL＝S1とする請求項１記載のデータプロセッサ。
【請求項４】
前記C1は前記データ処理ユニットの演算精度に対して丸めを必要としない正又は負の値である請求項２又は３記載のデータプロセッサ。
【請求項５】
前記関数は偶関数又は奇関数である請求項１記載のデータプロセッサ。
【請求項６】
コンピュータ装置が演算制御プログラムを実行することにより、所定の関数の値をその近似多項式F(x)を分解した和の多項式F1(x)+F2(x)を用いて取得するデータ処理方法であって、
前記式F1(x)は前記コンピュータ装置の演算精度に対して正確に計算可能な次数と係数を有しF1(x)＝BIG＋SMALLとして表現可能であり、F1(x)の絶対値はF2(x)の絶対値よりも大きく、BIGの絶対値はSMALLの絶対値よりも大きく、
前記多項式F1(x)＋F2(x)を用いた演算において、F(x)＝BIG＋SMALL＋F2(x)を求めるとき、BIG＋(SMALL＋F2(ｘ))を行うデータ処理方法。
【請求項７】
前記所定の関数の近似多項式は、F(x)＝x＋c1*x^3＋c2*x^5＋c3*x^7＋c4*x^9＋c5*x^11であり、
これに対応する和の多項式はF(x)＝(x＋C1*x^3)＋((c1−C1)*x^3＋c2*x^5＋c3*d1^7＋c4*x^9＋c5*x^11)であり、
F1(x)＝(x＋C1*x^3)であり、
F2(x)＝((c1−C1)*x^3＋c2*x^5＋c3*d1^7＋c4*x^9＋c5*x^11)であり、
C1は前記演算精度の範囲内のビット数で完全に表現できる浮動小数点数であり、
F1(x)/x＝(1＋C1*x^2)の部分演算において、1＋C1*x^2の演算結果の上位桁の値をB1に格納し、下位桁の値をS1に格納し、
BIG＝x*B1、SMALL＝x*S1とする請求項５記載のデータ処理方法。
【請求項８】
前記所定の関数の近似多項式は、F(x)＝c0＋c1*x^2＋c2*x^4＋c3*x^6＋c4*x^8＋c5*x^10であり、
これに対応する和の多項式はF(x)＝(c0＋C1*x^2)＋((c1−C1)*x^2＋c2*x^4＋c3*d1^6＋c4*x^8＋c5*x^10)であり、
F1(x)＝(c0＋C1*x^2)であり、
F2(x)＝((c1−C1)*x^2＋c2*x^4＋c3*d1^6＋c4*x^8＋c5*x^10)であり、
C1は前記演算精度の範囲内のビット数で完全に表現できる浮動小数点数であり、
F1(x)＝(c0＋C1*x^2)の部分演算において、c0＋C1*x^2の演算結果の上位桁の値をB1に格納し、下位桁の値をS1に格納し、
BIG＝B1、SMALL＝S1とする請求項５記載のデータ処理方法。
【請求項９】
コンピュータ装置によって実行される演算制御プログラムであって、前記演算制御プログラムは、所定の関数の値をその近似多項式F(x)を分解した和の多項式F1(x)+F2(x)を用いて取得するデータ処理の制御記述を有し、
前記データ処理に用いる前記式F1(x)は前記データ処理ユニットの演算精度に対して正確に計算可能な次数と係数を有しF1(x)＝BIG＋SMALLとして表現可能であり、F1(x)の絶対値はF2(x)の絶対値よりも大きく、BIGの絶対値はSMALLの絶対値よりも大きく、
前記データ処理は、前記多項式F1(x)＋F2(x)を用いた演算において、F(x)＝BIG＋SMALL＋F2(x)を求めるとき、GIB＋(SMALL＋F2(ｘ))を行う処理を含む演算制御プログラム。
【請求項１０】
前記所定の関数の近似多項式は、F(x)＝x＋c1*x^3＋c2*x^5＋c3*x^7＋c4*x^9＋c5*x^11であり、
これに対応する和の多項式はF(x)＝(x＋C1*x^3)＋((c1−C1)*x^3＋c2*x^5＋c3*d1^7＋c4*x^9＋c5*x^11)であり、
F1(x)＝(x＋C1*x^3)であり、
F2(x)＝((c1−C1)*x^3＋c2*x^5＋c3*d1^7＋c4*x^9＋c5*x^11)であり、
C1は前記演算精度の範囲内のビット数で完全に表現できる浮動小数点数であり、
F1(x)/x＝(1＋C1*x^2)の部分演算において、1＋C1*x^2の演算結果の上位桁の値をB1に格納し、下位桁の値をS1に格納し、
BIG＝x*B1、SMALL＝x*S1とする請求項９記載の演算制御プログラム。
【請求項１１】
前記所定の関数の近似多項式は、F(x)＝c0＋c1*x^2＋c2*x^4＋c3*x^6＋c4*x^8＋c5*x^10であり、
これに対応する和の多項式はF(x)＝(c0＋C1*x^2)＋((c1−C1)*x^2＋c2*x^4＋c3*d1^6＋c4*x^8＋c5*x^10)であり、
F1(x)＝(c0＋C1*x^2)であり、
F2(x)＝((c1−C1)*x^2＋c2*x^4＋c3*d1^6＋c4*x^8＋c5*x^10)であり、
C1は前記演算精度の範囲内のビット数で完全に表現できる浮動小数点数であり、
F1(x)＝(c0＋C1*x^2)の部分演算において、c0＋C1*x^2の演算結果の上位桁の値をB1に格納し、下位桁の値をS1に格納し、
BIG＝B1、SMALL＝S1とする請求項９記載の演算制御プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、sin(x)、atan(x)、sqrt(x)などの関数の値を近似多項式を用いて演算する技術に係り、例えばそのような関数の演算命令を実行可能なマイクロコンピュータに適用して有効な技術に関する。
【背景技術】
【０００２】
奇関数の近似多項式は例えば、F(x)＝c0*x＋c1*x^3＋c2*x^5＋c3*x^7（*は乗算記号、^は累乗記号である）となる。sin(x)、tan(x)などのように原点位置の微分係数が１となる関数の場合、最も低次のxの係数c0が１若しくは極端に１に近い値（単精度浮動小数点数では1.0fと記す）になる。このような近似多項式をホーナー法で計算する場合には、前記多項式近似をF(x)＝x(1＋x^2（c1＋x^2(c2＋x^2(c3))）)として、変数の高次側より積和演算を繰り返す。例えば最初にテンポラリな変数ansにc3を代入し（ans＝c3）、次にansにc2＋x^2*ansを代入する（ans＝c2＋x^2*ans）というように計算を行う。最後の積和演算部分では、ans＝1.0f＋ans*d2…(step1)、ans＝ans*d1…(step2)、が行われる。前記d2はd1*d1であり、d1は入力xの値である。
【０００３】
なお、多項式近似について記載された文献として非特許文献１、２がある。またＣ言語について記載された文献として非特許文献３がある。
【０００４】
【非特許文献１】CODY, W. J., AND WAITE, W. Software Manual for the Elementary Functions. Prentice-Hall, Englewood Cliffs, N.J., 1980.
【非特許文献２】Fike, C. T. 1968, Computer Evaluation of Mathematical Function(Englewood Cliffs, Nj: Prentice-Hall)
【非特許文献３】ISO/IEC: International Standard ISO/IEC 9899, Programming LanguageC, 1999
【発明の開示】
【発明が解決しようとする課題】
【０００５】
本発明者は関数の近似多項式による演算誤差について検討した。計算機で扱う実数に近い数値は浮動小数点数である。これは数学的な意味の実数とはまったく別物である。例えば単精度浮動小数点数(以下floatと呼ぶ)は32ビットのレジスタなどの入れ物に格納される。従って区間［−∞,∞］の値の内の2^32個の値だけが正確に表現できることになる。その他の値は近似値で我慢する他はない。例えば近似多項式F(x)＝x＋1/3.0f*x^3＋…について考える。初項xは入力であり誤差を持たない。しかし２番目の項は必ず誤差を持っている。まず、float定数1/3.0fは実数の1/3とは異なっている。この時点で誤差を持つ。さらに1/3.0f*x^3で誤差が拡大する。以下同様。しかし、入力xの絶対値が小さい間は、２番目の項の誤差は最終結果に殆ど影響を与えない。誤差を持たない初項xが２番目の項よりも格段に大きくなり、演算後の丸めで丸め方向を間違えることは稀にしか起こらないからである。入力xの絶対値が１に近くなれば、２番目の項の誤差は最終結果にかなり影響を与える。演算後の丸めで丸め方向をかなりの割合で間違えさせることになる。さらに入力xの絶対値が１を超えると、２番目の項の誤差が最終結果にそのまま反映されだす。この場合は最後の演算の桁落ちで、演算結果の精度(相対誤差の意味)が殆ど保証できなくなる。
【０００６】
本発明の目的は、演算装置による有効な演算ビット数を拡張することなく関数の近似多項式による演算精度を向上(丸め方向間違いの頻度を減少)させることにある。
【０００７】
本発明の前記並びにその他の目的と新規な特徴は本明細書の記述及び添付図面から明らかになるであろう。
【課題を解決するための手段】
【０００８】
本願において開示される発明のうち代表的なものの概要を簡単に説明すれば下記の通りである。
【０００９】
〔１〕データプロセッサはデータ処理ユニットと、前記データ処理ユニットが実行する演算制御プログラムを保有するプログラムメモリと、前記演算制御プログラムにしたがって前記データ処理ユニットがアクセスするワークメモリと、を有する。前記データ処理ユニットは、前記演算制御プログラムに従って所定の関数の値を演算するとき、近似多項式F(x)を分解した和の多項式F1(x)+F2(x)を用いる。前記式F1(x)は前記データ処理ユニットの演算精度に対して正確に計算可能な次数と係数を有しF1(x)＝BIG＋SMALLとして表現可能である。F1(x)の絶対値はF2(x)の絶対値よりも大きい。BIGの絶対値はSMALLの絶対値よりも大きい。前記多項式F1(x)＋F2(x)を用いた演算において、F(x)＝BIG＋SMALL＋F2(x)を求めるとき、BIG＋(SMALL＋F2(ｘ))の演算を行う。これにより、絶対値の小さいもの同士の加算で生じた丸め誤差はBIGに対して1度だけ作用する。丸め誤差が累積して大きな演算誤差を生ずる方向に最終解が丸められる可能性を低減することができる。
【００１０】
本発明の一つの具体的な形態として、前記所定の関数の近似多項式を、F(x)＝x＋c1*x^3＋c2*x^5＋c3*x^7＋c4*x^9＋c5*x^11とする。これに対応する和の多項式を、F(x)＝(x＋C1*x^3)＋((c1−C1)*x^3＋c2*x^5＋c3*d1^7＋c4*x^9＋c5*x^11)とする。F1(x)＝(x＋C1*x^3)であり、F2(x)＝((c1−C1)*x^3＋c2*x^5＋c3*d1^7＋c4*x^9＋c5*x^11)である。C1は前記演算精度の範囲内のビット数で完全に表現できる浮動小数点数である。F1(x)/x＝(1＋C1*x^2)の部分演算において、1＋C1*x^2の演算結果の上位桁の値をB1に格納し、下位桁の値をS1に格納する。S1の値はB1=1+C1*x^2の演算で失われた（落とされた）下位側ビットによる値であり、これによってそのような失われた値を回復している。そして、BIG＝x*B1、SMALL＝x*S1とする。これにより、BIGの絶対値はSMALLの絶対値よりも大きいこと、そして、F1(x)は前記データ処理ユニットの演算精度に対して正確に計算可能であることを保証することができる。
【００１１】
本発明の別の一つの具体的な形態として、前記所定の関数の近似多項式を、F(x)＝c0＋c1*x^2＋c2*x^4＋c3*x^6＋c4*x^8＋c5*x^10とする。これに対応する和の多項式をF(x)＝(c0＋C1*x^2)＋((c1−C1)*x^2＋c2*x^4＋c3*d1^6＋c4*x^8＋c5*x^10)とする。F1(x)＝(c0＋C1*x^2)であり、F2(x)＝((c1−C1)*x^2＋c2*x^4＋c3*d1^6＋c4*x^8＋c5*x^10)である。C1は前記演算精度の範囲内のビット数で完全に表現できる浮動小数点数である。F1(x)＝(c0＋C1*x^2)の部分演算において、1＋C1*x^2の演算結果の上位桁の値をB1に格納し、下位桁の値をS1に格納する。S1の値はB1=1+C1*x^2の演算で失われた（落とされた）下位側ビットによる値であり、これによってそのような失われた値を回復している。そして、BIG＝B1、SMALL＝S1とする。これにより、BIGの絶対値はSMALLの絶対値よりも大きいこと、そして、F1(x)は前記データ処理ユニットの演算精度に対して正確に計算可能であることを保証することができる。
【００１２】
〔２〕データ処理方法は、コンピュータ装置が演算制御プログラムを実行することにより、所定の関数の値をその近似多項式F(x)を分解した和の多項式F1(x)+F2(x)を用いて取得する。前記式F1(x)は前記コンピュータ装置の演算精度に対して正確に計算可能な次数と係数を有しF1(x)＝BIG＋SMALLとして表現可能である、F1(x)の絶対値はF2(x)の絶対値よりも大きく、BIGの絶対値はSMALLの絶対値よりも大きい。前記多項式F1(x)＋F2(x)を用いた演算において、F(x)＝BIG＋SMALL＋F2(x)を求めるとき、BIG＋(SMALL＋F2(ｘ))を行う。これにより、絶対値の小さいもの同士の加算で生じた丸め誤差はBIGに対して1度だけ作用する。丸め誤差が累積して大きな演算誤差を生ずる方向に最終解が丸められる可能性を低減することができる。
【００１３】
〔３〕コンピュータ装置によって実行される演算制御プログラムは、所定の関数の値をその近似多項式F(x)を分解した和の多項式F1(x)+F2(x)を用いて取得するデータ処理の制御記述を有する。前記データ処理に用いる前記式F1(x)は前記データ処理ユニットの演算精度に対して正確に計算可能な次数と係数を有しF1(x)＝BIG＋SMALLとして表現可能である。F1(x)の絶対値はF2(x)の絶対値よりも大きく、BIGの絶対値はSMALLの絶対値よりも大きい。前記データ処理は、前記多項式F1(x)＋F2(x)を用いた演算において、F(x)＝BIG＋SMALL＋F2(x)を求めるとき、GIB＋(SMALL＋F2(ｘ))を行う処理を含む。これにより、絶対値の小さいもの同士の加算で生じた丸め誤差はBIGに対して1度だけ作用する。丸め誤差が累積して大きな演算誤差を生ずる方向に最終解が丸められる可能性を低減することができる。
【発明の効果】
【００１４】
本願において開示される発明のうち代表的なものによって得られる効果を簡単に説明すれば下記の通りである。
【００１５】
すなわち、演算装置による有効な演算ビット数を拡張することなく関数の近似多項式による演算精度を向上させることができる。
【発明を実施するための最良の形態】
【００１６】
図１には本発明の一例に係るマイクロコンピュータが示される。同図に示されるマイクロコンピュータ１は、特に制限されないが、公知の半導体集積回路製造技術によって単結晶シリコンのような１個の半導体基板（半導体チップ）に形成される。このマイクロコンピュータ１は、データ処理ユニット２０として、例えば、整数演算を行なう中央処理装置（ＣＰＵ）２と共に、浮動小数点演算を行なう浮動小数点演算装置（ＦＰＵ）３を有する。更に、ディジタル信号処理で多用される積和演算に特化した積和演算装置４を有する。
【００１７】
中央処理装置２は内部アドレスバス１４及び内部データバス１３に結合される。中央処理装置２は、特に制限されないが、汎用レジスタや算術論理演算器で代表される演算部２Ａを有する。更に、中央処理装置２は、プログラムカウンタなどの制御用レジスタ群、そして命令のフェッチや解読並びに命令実行手順を制御したり演算制御を行う命令制御部２Ｂなどを有する。前記内部バス１３，１４には、前記ＣＰＵ２が実行する演算制御プログラムなどを保有するプログラムメモリ（ＰＧＭＭ）２１と、前記演算制御プログラムにしたがって前記ＣＰＵ２がアクセスするワークメモリ（ＷＲＫＭ）２２とを有する。中央処理装置２は、内部バス１３，１４に接続されたプログラムメモリ２１から命令をフェッチし、その命令を解読し、解読結果に応ずる制御信号を生成することにより、当該命令に応じたデータ処理を行う。
【００１８】
浮動小数点演算装置（ＦＰＵ）３及び積和演算装置４は内部データバス１３に結合される。浮動小数点演算装置（ＦＰＵ）３及び積和演算装置４は図示を省略する演算回路と共にデータレジスタを有し、このデータレジスタにメモリから演算データがロードされ、演算結果データは、そのレジスタからメモリにストアされる。前記ロード、ストアなどのためのアドレッシング動作はＣＰＵ２が行なう。したがって、ＦＰＵ３及び積和演算装置４はメモリアクセスのためのメモリアドレシング能力を備える必要はない。これは、ＦＰＵ３及び積和演算装置４によるメモリアドレシング回路の必要性を取り除いてチップ面積を節約するためである。
【００１９】
内部データバス１３及び内部アドレスバス１４はバスステートコントローラ６に結合される。マイクロコンピュータ１による外部アクセスは、前記バスステートコントローラ６に接続された外部バスインタフェース回路８で行う。外部バスインタフェース回路８は外部データバス１８及び外部アドレスバス１７に接続される。また、前記バスステートコントローラ６には、周辺データバス１６及び周辺アドレスバス１５を介して、例えば、クロック発生回路７、システムコントローラ１２、シリアルコミュニケーションインタフェースコントローラ（ＳＣＩ）９、タイマ１０及びＡ／Ｄコンバータ１１が結合される。それら周辺回路はデータレジスタや制御レジスタを有し、それらレジスタは前記バスステートコントローラ６を介してＣＰＵ２によってアクセスされる。
【００２０】
ＣＰＵ２が管理する内部メモリ空間、外部メモリ空間、前記ＳＣＩ９，タイマ１０，Ａ／Ｄコンバータ１１などの周辺回路に対するアドレスエリアの割り当ては予め決定されている。前記バスステートコントローラ６は、アクセスエリア毎にアクセスサイクル数やバス幅などがＣＰＵ２によって設定される図示を省略するバスコントロールレジスタを有し、ＣＰＵ２からのアクセスアドレスで指定されるメモリエリアに対するバス幅やアクセスサイクル数などのバス制御を行なって、バスサイクルを起動する。ＣＰＵ２からのアクセスの指示はバスコマンド２３としてバスステートコントローラ６に与えられる。バスコマンド２３にはアクセスサイズの指定や、リード、ライト、メモリアクセス等のストローブ信号が含まれる。
【００２１】
割込みコントローラ５は複数の割込み要求に対する優先制御やマスク制御などを行なって割込み信号２５をＣＰＵ２に与える。割込み要求は、バスステートコントローラ６からのバスエラー信号２４のほか、ＳＣＩ９、タイマ１０、Ａ／Ｄコンバータ１１そして外部からの図示を省略する割込み要求信号によって与えられる。
【００２２】
前記マイクロコンピュータ１は、特に制限されないが、クロック発生回路７から出力されるクロック信号ＣＬＫ０〜ＣＬＫ２に同期動作される。前記システムコントローラ１２はマイクロコンピュータ１の内蔵モジュールに対する動作停止の制御などを行う。
【００２３】
前記データ処理ユニット２０は、プログラムメモリ２１が保有する前記演算制御プログラムに従って所定の関数の値を演算するとき、その近似多項式F(x)を分解した和の多項式F1(x)+F2(x)を用いる。この多項式近似による演算として、奇関数sin(x)に適用した例を説明する。以下の説明において浮動小数点数は単精度とする。
【００２４】
前記近似多項式は、
sin(x) ≒F(x) = x ＋c1*x^3＋c2*x^5 ＋c3*x^7＋c4*x^9＋c5*x^11
とする。この近似多項式の係数は図２の値とされる。ただしErrは最大相対誤差である。
【００２５】
前記近似多項式F(x)を２つの多項式に分解し、
sin(x)＝F(x)＝(x ＋C1*x^3) ＋((c1−C1)*x^3 ＋c2*x^5＋c3*x^7＋…)＝F1(x)＋F2(x)
とする。C1は−0.149335208038861fであり、単精度浮動小数点数によって誤差なく表現できる値とする。またF1(x)の絶対値はF2(x)の絶対値よりも大きく、例えばF2(x)/F(x)の絶対値は最大で0.01程度である。要するに、前記式F1(x)は前記データ処理ユニットの演算精度に対して正確に計算可能な次数と係数を有する。そして、F1(x)＝BIG＋SMALLとして表現する。BIGの絶対値はSMALLの絶対値よりも大きくする。前記多項式F1(x)＋F2(x)を用いた演算において、F(x)＝BIG＋SMALL＋F2(x)を求めるとき、GIB＋(SMALL＋F2(ｘ))の演算を行う。これにより、絶対値の小さいもの同士の加算で生じた丸め誤差はBIGに対して1度だけ作用する。丸め誤差が累積して大きな演算誤差を生ずる方向に最終解が丸められる可能性を低減することができる。
【００２６】
図３には上記多項式演算による演算方法の具体例がC言語記述で示される。Floatは単精度浮動小数点数を意味する。入力xをd1、テンポラリな変数をd2,d2f,d3,b1,s1とする。b1はBIG/x、s1はSMALL/xである。
【００２７】
図３において、step1,step2,step3で多項式F1(d1)/d1=1.0+C1*d1^2を精度良く計算する。即ち、step1でd1^2を精度良く計算する。d2=d1*d1の計算で失われた小さい値をd2f=d1*d1-d2で回復している。この演算は前記非特許文献３における積和演算fmaf()を前提にしており、積の値総てが加算(減算)に参加することになる。
【００２８】
Step2ではb1=1.0f+C1*d2を精度良く計算する。b1=1.0f＋C1*d2で失われた小さい値をs1=1.0f−b1、s1−=C1*d2の２つの実行文で回復している。
【００２９】
1.0f+C1*d2== b1 +s1
ここでの積和演算も前記非特許文献３における積和演算fmaf()を前提にしている。Step3でd1^2とd2との差分であるd2fの寄与をs1に加算している。
【００３０】
1.0f+C1*d1*d1==1.0f+C1*(d2+d2f)==1.0f+C1*d2+C1*d2f==b1+(s1+C1*d2f)=b1+新s1
以上でF1(d1)/d1=1.0+C1*d1*d1を精度良く演算したことになる。要するに、1.0+C1*d1*d1 ＝＝ b1 + s1ということになる。但し、単精度浮動小数点数を用いているので、上記右辺のb1 + s1をこの時点で計算すればs1の情報が失われてしまうので、この段階では右辺の演算b1 + s1を留保しておく必要がある。
【００３１】
図３のstep4からstep5まででF2(x)/xを計算している。ここではホーナー法による計算を行っている。Step6からstep7の演算ではsin(x)の近似多項式F(x)を以下の式で表されるようにして演算している。即ち、Sin(d1)≒F(d1)＝F1(d1)＋F2(d1)＝d1*( b1＋s1) ＋d1^3*ans…（ansはstep5直後のもの）とし、更にF(d1)＝d1*b1＋(d1*s1＋d1^3*ans)）とする。この最後の式通りに計算すれば、最終結果の丸め方向間違いを最小化できる。初項の絶対値が圧倒的に大きく、しかも誤差を持たないからである。最後の積和演算も前記非特許文献３における積和演算fmaf()を前提にしている。
【００３２】
図４にはsin(x)の近似多項式に対して単なるホーナー法を用いて演算を行う従来からの一般的な演算手順が示される。奇関数sin(x)の区間[-1.6, 1.6]を近似する多項式F(x) = x + c1*x^3 + c2*x^5+...+c5*x^11
をそのままホーナー法で計算する例である。
【００３３】
図５には３種類の異なる演算区間において図４の従来からの演算手法と本発明の演算手法とによって生ずる丸め方向誤りの発生率が対比される。丸め方向に誤りを生じているかを判定するためなお基準の値は、倍精度で計算したsin(x)の値を単精度に丸めたものである。従来法では、入力の絶対値が大きくなると丸め方向を頻繁に間違える。しかし本発明の算法では入力の絶対値が大きい場合でも丸め方向の間違いは少ないままである。
【００３４】
本発明は偶関数にも適用可能である。例えばcos(x)の例を簡単に説明する。最初の２つをマクローリン展開の係数1.0、−0.5に固定した近似多項式を導出できる。そして下記式
cos(x) ≒ F(x)＝1.0−0.5*x^2 ＋c2*x^4＋c3*x^6＋…＝F1(x)＋F2(x)
に示されるように、最初の２つの項をF1(x)とする。残りをF2(x)とする。このF1(x)の計算は図３のstep1、step2、step3で実現できる。このF2(x)の計算は図３のstep4からstep5と同様に計算できる。そして以下の等式
F(x) ＝b1＋s1＋F2(x)＝b1＋(s1＋F2(x))
の最後の式通りに計算すれば、丸め方向間違いを最小化することができる。
【００３５】
偶関数cos(x)の更に別の例を説明する。最初の２つをマクローリン展開の係数1.0、−0.5に固定した近似多項式を導出できる。そしてF1(x)を以下の式
F1(x)＝(1.0−0.25*x^2)^2
の通りとする。そして残りをF2(x)とする。要するに、下記式
cos(x)≒F(x)＝1.0−0.25*x^2)^2＋(c2-0.25^2)*x^4＋c3*x^6＋…＝ F1(x)＋F2(x)
のようにする。このF1(x)の計算は図３のstep1,step2,step3でF1(x)の平方根をb1,s1に求めておき、以下の近似式
F1(x)＝(b1＋s1)^2≒b1*b1＋2*b1*s1
で実現できる。なお、最初の項の乗算と２つの項の加算は留保する。そしてF2(x)の計算は図３のstep4からstep5と同様に計算できる。そして以下の等式
F(x)＝b1*b1＋2*b1*s1＋F2(x)＝b1*b1＋(2*b1*s1＋F2(x))
の最後の式通りに計算すれば、丸め方向間違いを最小化することができる。
【００３６】
以上本発明者によってなされた発明を実施形態に基づいて具体的に説明したが、本発明はそれに限定されるものではなく、その要旨を逸脱しない範囲において種々変更可能であることは言うまでもない。
【００３７】
例えば、例えば本発明に適用可能な関数はsin(x)、cos(x)に限定されず、atan(x)など種々の奇関数及び偶関数に関する近似多項式演算に適用することができる。またlog(1+x)にも適用できる。
【図面の簡単な説明】
【００３８】
【図１】本発明の一例に係るマイクロコンピュータのブロック図である。
【図２】sin(x)の近似多項式F(x) = x ＋c1*x^3＋c2*x^5 ＋c3*x^7＋c4*x^9＋c5*x^11における係数の説明図である。
【図３】本発明に係る記多項式演算による演算方法の具体例をC言語記述で示す説明図である。
【図４】sin(x)の近似多項式に対して単なるホーナー法を用いて演算を行う従来からの一般的な演算方法をC言語標記で示す説明図である。
【図５】３種類の異なる演算区間において図４の従来からの演算手法と図３の本発明の演算手法とによって生ずる丸め方向誤りの発生率を対比して示す説明図である。
【符号の説明】
【００３９】
１  マイクロコンピュータ
２  中央処理装置
３  浮動小数点演算装置
４  積和演算装置
２０  データ処理ユニット
２１  プログラムメモリ
２２  ワークメモリ
【図１】
2006323710.tif 000003
【図２】
2006323710.tif 000004
【図３】
2006323710.tif 000005
【図４】
2006323710.tif 000006
【図５】
2006323710.tif 000007