除算アルゴリズムについて

質問

1/x　を、乗算器しかないDSPなどで演算するアルゴリズムについて
お教え願えないでしょうか？
おぼろげな記憶では、0＜ｘ≦1　などのように正規化して、
バタフライ演算を行って求めるアルゴリズムがあったと記憶しますが、
正確な事がわかりません。
宜しくお願いします。

回答No.3

stomachman · ベストアンサー

No.1のコメントに対する回答です。

f[1]=f[0](2-xf[0]) = (1-e)(1+e)/x = (1-e^2)/x
というのはつまり、誤差eの１次の項をうち消すようにしてやればよいのです。「簡単に計算できる関数の逆関数」あるいは一般に「簡単に計算できる検算法がある関数」を計算するときに広く使えます。

この漸化式はニュートン法から導くこともできます。
g(f)=x-1/f
とおいて、方程式
g(f)=0
を解く。ニュートン法を使うと、
g'(f) = dg/df = 1/(f^2)
を用いて
f[n+1]=f[n]-g(f[n])/g'(f[n])
= f[n]-(x-1/f[n])/(1/(f[n]^2))
= f[n]-(x(f[n]^2)-f[n])
= 2f[n]-x(f[n]^2)
= f[n](2-xf[n])

　この場合は旨く行きましたが、ニュートン法がわり算なしの漸化式になるとは限りません。そういうときは、わり算なしで誤差|e|が繰り返しの度に単調減少するような計算を何らか工夫してやることになります。

　なお、実際の応用においてはtable lookupは馬鹿にならない性能を発揮します。10bitの精度で1/100をやるんだったら、高々数百個の要素を持つテーブルを用意すれば十分で、２次あるいは３次のラグランジュ補間を使えば要素をさらに１桁以上減らせます。メモリやキャッシュの隅っこに十分収まってしまいますね。

お礼

stomachmanさま、お世話になってます。
成る程～～　　反復法やニュートンラフソン法　の最初のところに
f[n+1]=f[n]-g(f[n])/g'(f[n])
が載っていますね（＾＾；）
g(f)=x-1/f 　を　fx-1　と変形してしまったために
代入すると全部消えてしまって導けなかった（？）ようです。

実のところ、分母には高々数十種類（6bitで充分）の整数と決まっているので
全部テーブル引きでも良いのです（実際、テーブルを小さくしてラグランジュ補間
を実装するのはかなり大変だろうと思います）が、もう少し方法はないものかと
思い、ここに投稿した次第です。
親切に解説していただき有難うございました。

また、自己レスにはなりますが、
DesignWaveマガジン　1999年12月号、2000年5月号に、全く同じアルゴリズムで
除算器を構成する方法が記載されていました。

回答No.2

cherry_moon

固定小数点演算なのか、浮動小数点演算なのか。
除算器を作ろうとしてるのか、ソフトウェアで計算したいのか。

とかの情報も書いた方がいいですよ。

浮動小数点演算をソフトウェアで実装するなら、
stomachman さんの方法でいいと思います。

お礼

cherry_moonさま、お世話になります（＾＾）
情報はなるべく多いほうが良いというわけですね。
ありがとうございます。

状況を説明しますと、
１　固定小数点で良い
２　ハードウェア（FPGA）で除算器もどきを実装したい
３　単純に引いた回数をカウントしたのでは間に合わない
４　大雑把にいうと10万/100　程度の演算です。
５　４の整数解が得られれば良い
６　実は精度も0.1％（10bit）程度で良い
なので、最も適した回路を模索しているところです。

stomachmanさんの説明のアルゴリズムから大雑把に規模・スピードを見積もりましたが、
結構小規模で速度・精度が得られそうです。

回答No.1

stomachman

1/xの近似値をf[0]とします。その相対誤差をe（|e|<1)とすると
f[0] = (1+e)/x
という関係になっています。
次に
f[1]=f[0](2-xf[0])
を計算する。
2-xf[0]= 2-(1+e) = (1-e)
ですから、
f[1]=f[0](2-xf[0]) = (1-e)(1+e)/x = (1-e^2)/x
従って誤差はeだったものがe^2に改良されました。これを
f[n+1]=f[n](2-xf[n])
のように繰り返せば、誤差はe→e^2→e^4→e^8…とどんどん小さくなり、有効桁数が繰り返しの度に倍になっていきます。

例えばx=0.1の逆数を求めるのに出発値として1/x ≒ f[0] =2-xを使ったとすると、

f[0]=1.9
f[1]=3.439
f[2]=5.6953279
f[3]=8.146979811
f[4]=9.656631618
f[5]=9.988209815
f[6]=9.999986099

どのぐらいの範囲のXが入ってくるかによって、前処理の方法がいろいろあり得ます。また必要な答の精度が荒くて良いのなら表を引く(table lookup)、表を引いてから２次補間する、などの方法でも十分ですし、表を出発値f[0]を求めるのに使うのも良い方法です。

お礼

はじめまして、stomachmanさま。ありがとうございます。
成る程、これは使えそうです。
あまえついでなのですが、
f[n+1]=f[n](2-xf[n])
特性方程式？　はどうやって導くのでしょうか？
初歩的な数値計算法の本をいくつかあさってみましたが
恥ずかしながら判りませんでした。