5

リファレンス実装として、次のことを考慮してください。

/* calculates (a * b) / c */
uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c)
{
    uint64_t x = a;
    x = x * b;
    x = x / c;
    return x;
}

64ビット整数型を必要としない実装（Cまたは擬似コード）に興味があります。

私は次のような概要を示す実装のスケッチを開始しました。

/* calculates (a * b) / c */
uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c)
{
    uint32_t d1, d2, d1d2;
    d1 = (1 << 10);
    d2 = (1 << 10);
    d1d2 = (1 << 20); /* d1 * d2 */
    return ((a / d1) * (b /d2)) / (c / d1d2);
}

ただし、オーバーフローを回避し（（a / d1）*（b / d2）<= UINT32_MAX）、計算全体のエラーを最小限に抑えることができるd1とd2の値を選択するのは困難です。

何かご意見は？

4

7 に答える 7

4

Paulによって投稿されたアルゴリズムをunsignedintに適合させました（符号を処理する部分を省略しています）。アルゴリズムは基本的に古代エジプトの分数との乗算です（スラッシュはここで実際の除算を示します）。afloor(b/c) + (b%c)/c

uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c)
{
    uint32_t q = 0;              // the quotient
    uint32_t r = 0;              // the remainder
    uint32_t qn = b / c;
    uint32_t rn = b % c;
    while(a)
    {
        if (a & 1)
        {
            q += qn;
            r += rn;
            if (r >= c)
            {
                q++;
                r -= c;
            }
        }
        a  >>= 1;
        qn <<= 1;
        rn <<= 1;
        if (rn >= c)
        {
            qn++; 
            rn -= c;
        }
    }
    return q;
}

このアルゴリズムは、32ビットに収まる限り、正確な答えを生成します。オプションで、余りを返すこともできますr。

于 2010-11-10T13:29:14.593 に答える
3

最も簡単な方法は、中間結果を64ビットに変換することですが、cの値に応じて、別のアプローチを使用できます。

((a/c)*b  +  (a%c)*(b/c) + ((a%c)*(b%c))/c

唯一の問題は、の値が大きい場合でも最後の項がオーバーフローする可能性があることですc。まだそれについて考えています。

于 2010-11-10T12:13:31.620 に答える
3

最初にaをcで除算し、除算のリマインダーを取得し、cで除算する前にリマインダーにbを掛けることができます。そうすれば、最後の除算でのみデータが失われ、64ビット除算を行った場合と同じ結果が得られます。

次のように式を書き直すことができます（ここで、\は整数除算です）。

a * b / c =
(a / c) * b =
(a \ c + (a % c) / c) * b =
(a \ c) * b + ((a % c) * b) / c

a> = bであることを確認することにより、オーバーフローする前に、より大きな値を使用できます。

uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c) {
  uint32_t hi = a > b ? a : b;
  uint32_t lo = a > b ? b : a;
  return (hi / c) * lo + (hi % c) * lo / c;
}

もう1つのアプローチは、乗算と除算の代わりに加算と減算をループすることですが、もちろんそれははるかに多くの作業です。

uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c) {
  uint32_t hi = a > b ? a : b;
  uint32_t lo = a > b ? b : a;
  uint32_t sum = 0;
  uint32_t cnt = 0;
  for (uint32_t i = 0; i < hi; i++) {
    sum += lo;
    while (sum >= c) {
      sum -= c;
      cnt++;
    }
  }
  return cnt;
}
于 2010-11-10T12:37:58.887 に答える
3

www.google.com/codesearchで検索すると、この素晴らしく明白なものを含む、多くの実装が見つかります。私は特に広範なコメントとよく選ばれた変数名が好きです

INT32 muldiv(INT32 a, INT32 b, INT32 c)
{ INT32 q=0, r=0, qn, rn;
  int qneg=0, rneg=0;
  if (c==0) c=1;
  if (a<0) { qneg=!qneg; rneg=!rneg; a = -a; }
  if (b<0) { qneg=!qneg; rneg=!rneg; b = -b; }
  if (c<0) { qneg=!qneg;             c = -c; }

  qn = b / c;
  rn = b % c;

  while(a)
  { if (a&1) { q += qn;
               r += rn;
               if(r>=c) { q++; r -= c; }
             }
    a  >>= 1;
    qn <<= 1;
    rn <<= 1;
    if (rn>=c) {qn++; rn -= c; }
  }
  result2 = rneg ? -r : r;
  return qneg ? -q : q;
}

http://www.google.com/codesearch/p?hl=en#HTrPUplLEaU/users/mr/MCPL/mcpl.tgz|gIE-sNMlwIs/MCPL/mintcode/sysc/mintsys.c&q=muldiv%20lang:c

于 2010-11-10T12:53:11.133 に答える
1

SvenのコードをUINT16として実装し、集中的にテストしました。

uint16_t muldiv16(uint16_t a, uint16_t b, uint16_t c);

int main(int argc, char *argv[]){
    uint32_t a;
    uint32_t b;
    uint32_t c;
    uint16_t r1, r2;

// ~167 days, estimated on i7 6700k, single thread.
// Split the 'a' range, to run several instances of this code on multi-cores processor
// ~1s, with an UINT8 implementation
    for(a=0; a<=UINT16_MAX; a++){
        for(b=0; b<=UINT16_MAX; b++){
            for(c=1; c<=UINT16_MAX; c++){
                r1 = uint16_t( a*b/c );
                r2 = muldiv16(uint16_t(a), uint16_t(b), uint16_t(c));
                if( r1 != r2 ){
                    std::cout << "Err: " << a << " * " << b << " / " << c << ", result: " << r2 << ", exected: " << r1 << std::endl;
                    return -1;
                }
            }
        }
        std::cout << a << std::endl
    }
    std::cout << "Done." << std::endl;
    return 0;
}

残念ながら、「b」（0-2147483647）はUINT31に制限されているようです。

これが私の修正です。これは機能しているようです（UINT16でのテストは完了していませんが、たくさん実行しています。UINT8で完了しました）。

uint32_t muldiv32(uint32_t a, uint32_t b, uint32_t c)
{
    uint32_t q = 0;              // the quotient
    uint32_t r = 0;              // the remainder
    uint32_t qn = b / c;
    uint32_t rn = b % c;
    uint32_t r_carry;
    uint32_t rn_carry;
    while(a)
    {
        if (a & 1)
        {
            q += qn;
            r_carry = (r > UINT32_MAX-rn);
            r += rn;
            if (r >= c || r_carry)
            {
                q++;
                r -= c;
            }
        }
        a  >>= 1;
        qn <<= 1;
        rn_carry = rn & 0x80000000UL;
        rn <<= 1;
        if (rn >= c || rn_carry)
        {
            qn++;
            rn -= c;
        }
    }
    return q;
}

編集：残りを返し、ラウンドを管理し、オーバーフローについて警告し、もちろん、a、b、およびcのUINT32の全範囲を管理する改善：

typedef enum{
    ROUND_DOWNWARD=0,
    ROUND_TONEAREST,
    ROUND_UPWARD
}ROUND;

//remainder is always positive for ROUND_DOWN ( a * b = c * q + remainder )
//remainder is always negative for ROUND_UPWARD ( a * b = c * q - remainder )
//remainder is signed for ROUND_CLOSEST ( a * b = c * q + sint32_t(remainder) )
uint32_t muldiv32(uint32_t a, uint32_t b, uint32_t c, uint32_t *remainder, ROUND round, uint8_t *ovf)
{
    uint32_t q = 0;              // the quotient
    uint32_t r = 0;              // the remainder
    uint32_t qn = b / c;
    uint32_t rn = b % c;
    uint32_t r_carry;
    uint32_t rn_carry;
    uint8_t o = 0;
    uint8_t rup;
    while(a)
    {
        if (a & 1)
        {
            o |= (q > UINT32_MAX-qn);
            q += qn;
            r_carry = (r > UINT32_MAX-rn);
            r += rn;
            if (r >= c || r_carry)
            {
                o |= (q == UINT32_MAX);
                q++;
                r -= c;
            }
        }
        a  >>= 1;
        qn <<= 1;
        rn_carry = rn & 0x80000000;
        rn <<= 1;
        if (rn >= c || rn_carry)
        {
            qn++;
            rn -= c;
        }
    }
    rup = (round == ROUND_UPWARD && r);
    rup |= (round == ROUND_TONEAREST && ((r<<1) >= c || r & 0x80000000));
    if(rup)
    {   //round
        o |= (q == UINT32_MAX);
        q++;
        r = (round == ROUND_UPWARD) ? c-r : r-c;
    }
    if(remainder)
        *remainder = r;
    if(ovf)
        *ovf = o;
    return q;
}

おそらく、別のアプローチが存在する可能性があります。おそらくさらに効率的です。8ビット、16ビット、および32ビットのMCUは、64ビットの計算（long long int）を計算できます。コンパイラがそれをどのようにエミュレートするか知っている人はいますか？

編集2：

8ビットMCUでのいくつかの興味深いタイミングは次のとおりです。

UINT8 x UINT8 / UINT8：3.5µs

UINT16 x UINT16 / UINT16：22.5µs、muldiv8：29.9〜45.3µs

UINT32 x UINT32 / UINT32：84µs、muldiv16：120〜189µs

FLOAT32 * FLOAT32 / FLOAT32：40.2 ot 135.5µs、muldiv32：1.193〜1.764ms

そして32ビットMCUの場合：

タイプ-最適化されたコード-最適化なし

UINT32：521ns-604ns

UINT64：2958ns-3313ns

FLOAT32：2563ns〜2688ns

muldiv32：6791ns-25375ns

したがって、コンパイラはこのCアルゴリズムよりも賢いです。また、ネイティブレジスタよりも整数が大きい場合よりも（FPUがなくても）float変数を使用する方が常に優れています（float32の精度はuint32よりも低く、16777217以降です）。

Edit3：わかりました。私のNビットMCUはN-bits MUL N-bits、2Nビットの結果を生成するネイティブ命令を使用しており、2つのNビットレジスタに格納されています。

ここで、Cの実装を見つけることができます（EasyasPiのソリューションをお勧めします）

2N-bits DIV N-bitsしかし、彼らはネイティブの指示を持っていません。代わりに、ループと2Nビット変数（ここではUINT64）を使用して、gccの__udivdi3関数を使用しています。したがって、これは元の質問の解決策にはなりません。

于 2021-02-28T08:58:39.367 に答える
0

bとcが両方とも定数である場合、エジプト式分数を使用して非常に簡単に結果を計算できます。

例えば。y = a*4/99は次のように書くことができます

y = a / 25 + a / 2475

Cのエジプト式分数への回答で説明されているように、任意の分数をエジプト式分数の合計として表すことができます。

bとcを事前に修正しておくことは少し制限のように思えるかもしれませんが、この方法は他の人が答える一般的なケースよりもはるかに簡単です。

于 2015-11-18T11:35:38.483 に答える
-4

できない理由があると思います

x = a/c;
x = x*b;

ある？そして多分追加

y = b/c;
y = y*a;

if ( x != y )
    return ERROR_VALUE;

整数除算を使用しているためa*b/c、またはよりも大きいa/c*b場合は異なる値になる可能性があることに注意してください。また、との両方がそれよりも小さい場合は機能しません。cab abc

于 2010-11-10T12:34:15.170 に答える