Умножение чисел произвольной длинны взял из книги "Алгоритмические трюки для программистов" (Уоррен). Тут пришлось поковыряться немного, в моем скане куча ошибок в исходниках, да и стиль у автора не очень внятный.
#include <limits>
#include <algorithm>
inline bool unsigned_addition( UINT16 a, UINT16 b, UINT16 *c = 0 )
{
UINT16 r = a + b; if (c) *c = r; return !(r>=a && r>=b);
}
inline bool unsigned_addition( UINT32 a, UINT32 b, UINT32 *c = 0 )
{
UINT32 r = a + b; if (c) *c = r; return !(r>=a && r>=b);
}
inline bool unsigned_addition( UINT64 a, UINT64 b, UINT64 *c = 0 )
{
UINT64 r = a + b; if (c) *c = r; return !(r>=a && r>=b);
}
inline bool unsigned_subtraction( UINT16 a, UINT16 b, UINT16 *c = 0 )
{
UINT16 r = a - b; if (c) *c = r; return a<b;
}
inline bool unsigned_subtraction( UINT32 a, UINT32 b, UINT32 *c = 0 )
{
UINT32 r = a - b; if (c) *c = r; return a<b;
}
inline bool unsigned_subtraction( UINT64 a, UINT64 b, UINT64 *c = 0 )
{
UINT64 r = a - b; if (c) *c = r; return a<b;
}
inline bool unsigned_addition( UINT16 a, UINT16 b, UINT16 &c )
{ return unsigned_addition( a, b, &c ); }
inline bool unsigned_addition( UINT32 a, UINT32 b, UINT32 &c )
{ return unsigned_addition( a, b, &c ); }
inline bool unsigned_addition( UINT64 a, UINT64 b, UINT64 &c )
{ return unsigned_addition( a, b, &c ); }
inline bool unsigned_subtraction( UINT16 a, UINT16 b, UINT16 &c )
{ return unsigned_subtraction( a, b, &c ); }
inline bool unsigned_subtraction( UINT32 a, UINT32 b, UINT32 &c )
{ return unsigned_subtraction( a, b, &c ); }
inline bool unsigned_subtraction( UINT64 a, UINT64 b, UINT64 &c )
{ return unsigned_subtraction( a, b, &c ); }
// returns true if carry detectedinline bool signed_addition( INT16 a, INT16 b, INT16 *c = 0 )
{
INT16 r = a + b; if (c) *c = r; return !(b >= 0 == r >= a);
}
inline bool signed_addition( INT32 a, INT32 b, INT32 *c = 0 )
{
INT32 r = a + b; if (c) *c = r; return !(b >= 0 == r >= a);
}
inline bool signed_addition( INT64 a, INT64 b, INT64 *c = 0 )
{
INT64 r = a + b; if (c) *c = r; return !(b >= 0 == r >= a);
}
inline bool signed_subtraction( INT16 a, INT16 b, INT16 *c = 0 )
{ return signed_addition( a, -b, c ); }
inline bool signed_subtraction( INT32 a, INT32 b, INT32 *c = 0 )
{ return signed_addition( a, -b, c ); }
inline bool signed_subtraction( INT64 a, INT64 b, INT64 *c = 0 )
{ return signed_addition( a, -b, c ); }
inline bool signed_addition( INT16 a, INT16 b, INT16 &c )
{ return signed_addition( a, b, &c ); }
inline bool signed_addition( INT32 a, INT32 b, INT32 &c )
{ return signed_addition( a, b, &c ); }
inline bool signed_addition( INT64 a, INT64 b, INT64 &c )
{ return signed_addition( a, b, &c ); }
inline bool signed_subtraction( INT16 a, INT16 b, INT16 &c )
{ return signed_subtraction( a, b, &c ); }
inline bool signed_subtraction( INT32 a, INT32 b, INT32 &c )
{ return signed_subtraction( a, b, &c ); }
inline bool signed_subtraction( INT64 a, INT64 b, INT64 &c )
{ return signed_subtraction( a, b, &c ); }
template<typename IntT>
bool integer_addition( IntT a, IntT b, IntT *c = 0)
{
return signed_addition( (INT64)a, (INT64)b, (INT64*)c );
}
template<typename IntT>
bool integer_addition( IntT a, IntT b, IntT &c )
{
INT64 tmp;
bool bRes = signed_addition( (INT64)a, (INT64)b, &tmp );
c = (IntT)tmp;
return bRes;
}
template<> bool integer_addition<INT16> (INT16 a, INT16 b, INT16 *c) { return signed_addition( a, b, c ); }
template<> bool integer_addition<INT16> (INT16 a, INT16 b, INT16 &c) { return signed_addition( a, b, c ); }
template<> bool integer_addition<INT32> (INT32 a, INT32 b, INT32 *c) { return signed_addition( a, b, c ); }
template<> bool integer_addition<INT32> (INT32 a, INT32 b, INT32 &c) { return signed_addition( a, b, c ); }
template<> bool integer_addition<INT64> (INT64 a, INT64 b, INT64 *c) { return signed_addition( a, b, c ); }
template<> bool integer_addition<INT64> (INT64 a, INT64 b, INT64 &c) { return signed_addition( a, b, c ); }
template<> bool integer_addition<UINT16>(UINT16 a, UINT16 b, UINT16 *c) { return unsigned_addition( a, b, c ); }
template<> bool integer_addition<UINT16>(UINT16 a, UINT16 b, UINT16 &c) { return unsigned_addition( a, b, c ); }
template<> bool integer_addition<UINT32>(UINT32 a, UINT32 b, UINT32 *c) { return unsigned_addition( a, b, c ); }
template<> bool integer_addition<UINT32>(UINT32 a, UINT32 b, UINT32 &c) { return unsigned_addition( a, b, c ); }
template<> bool integer_addition<UINT64>(UINT64 a, UINT64 b, UINT64 *c) { return unsigned_addition( a, b, c ); }
template<> bool integer_addition<UINT64>(UINT64 a, UINT64 b, UINT64 &c) { return unsigned_addition( a, b, c ); }
template<typename IntT>
bool integer_subtraction( IntT a, IntT b, IntT *c = 0)
{
return signed_subtraction( (INT64)a, (INT64)b, (INT64*)c );
}
template<typename IntT>
bool integer_subtraction( IntT a, IntT b, IntT &c )
{
INT64 tmp;
bool bRes = signed_subtraction( (INT64)a, (INT64)b, &tmp );
c = (IntT)tmp;
return bRes;
}
template<> bool integer_subtraction<INT16> (INT16 a, INT16 b, INT16 *c) { return signed_subtraction( a, b, c ); }
template<> bool integer_subtraction<INT16> (INT16 a, INT16 b, INT16 &c) { return signed_subtraction( a, b, c ); }
template<> bool integer_subtraction<INT32> (INT32 a, INT32 b, INT32 *c) { return signed_subtraction( a, b, c ); }
template<> bool integer_subtraction<INT32> (INT32 a, INT32 b, INT32 &c) { return signed_subtraction( a, b, c ); }
template<> bool integer_subtraction<INT64> (INT64 a, INT64 b, INT64 *c) { return signed_subtraction( a, b, c ); }
template<> bool integer_subtraction<INT64> (INT64 a, INT64 b, INT64 &c) { return signed_subtraction( a, b, c ); }
template<> bool integer_subtraction<UINT16>(UINT16 a, UINT16 b, UINT16 *c) { return unsigned_subtraction( a, b, c ); }
template<> bool integer_subtraction<UINT16>(UINT16 a, UINT16 b, UINT16 &c) { return unsigned_subtraction( a, b, c ); }
template<> bool integer_subtraction<UINT32>(UINT32 a, UINT32 b, UINT32 *c) { return unsigned_subtraction( a, b, c ); }
template<> bool integer_subtraction<UINT32>(UINT32 a, UINT32 b, UINT32 &c) { return unsigned_subtraction( a, b, c ); }
template<> bool integer_subtraction<UINT64>(UINT64 a, UINT64 b, UINT64 *c) { return unsigned_subtraction( a, b, c ); }
template<> bool integer_subtraction<UINT64>(UINT64 a, UINT64 b, UINT64 &c) { return unsigned_subtraction( a, b, c ); }
template<typename IntT>
IntT integer_abs( IntT i ) { return (i<0) ? -i : return i; }
template<typename IntT>
IntT integer_sign( IntT i ) { return (i<0) ? -1 : 1; }
template<typename IntT>
int sign_multiple( IntT a, IntT b )
{ return integer_sign(a)==integer_sign(b) ? 1 : -1; }
template<typename IntT>
IntT make_lo_half_mask()
{
IntT mask = 0xFF;
IntT resMask = mask;
for(unsigned i =0; i!=(unsigned)(sizeof(IntT)/2); ++i)
{
resMask |= mask;
mask <<= 8;
}
return resMask;
}
template<typename IntT>
unsigned calc_half_shift() { return (unsigned)(sizeof(IntT)*4); }
template<typename IntT>
IntT lo_part( IntT i ) { return i & make_lo_half_mask<IntT>(); }
template<typename IntT>
IntT hi_part( IntT i )
{ return (i >> calc_half_shift<IntT>()) & make_lo_half_mask<IntT>(); }
/* keep sign */template<typename IntT>
IntT lo_part_s( IntT i ) { return i & make_lo_half_mask<IntT>(); }
template<typename IntT>
IntT hi_part_s( IntT i )
{ return i >> calc_half_shift<IntT>(); }
template<typename UIntT, typename UIntT2>
void multiply_num_mn_unsigned_impl( UIntT *w
, const UIntT *u, const UIntT *v
, unsigned m, unsigned n
)
{
UIntT2 k, t; //, b;unsigned i, j; // countersfor (i=0; i<m; i++)
w[i] = 0;
for(j=0; j!=n; j++)
{
k = 0;
for(i=0; i!=m; i++)
{
t = u[i]*v[j] + w[i+j] + k;
w[i+j] = static_cast<UIntT>(lo_part_s(t)); // t; // (Т.е., t & OxFFFF).
k = hi_part_s(t); // t >> 16;
}
w[j+m] = static_cast<UIntT>(k);
}
}
inline void
multiply_num_mn_unsigned( UINT8 *w, const UINT8 *u, const UINT8 *v
, unsigned m, unsigned n )
{ multiply_num_mn_unsigned_impl<UINT8, UINT16>(w, u, v, m, n ); }
inline void
multiply_num_mn_unsigned( UINT16 *w, const UINT16 *u, const UINT16 *v
, unsigned m, unsigned n )
{ multiply_num_mn_unsigned_impl<UINT16, UINT32>(w, u, v, m, n ); }
inline void
multiply_num_mn_unsigned( UINT32 *w, const UINT32 *u, const UINT32 *v
, unsigned m, unsigned n )
{ multiply_num_mn_unsigned_impl<UINT32, UINT64>(w, u, v, m, n ); }
template<typename UIntT, typename SIntT, typename UIntT2>
void multiply_num_mn_signed_impl( UIntT *w
, const UIntT *u, const UIntT *v
, unsigned m, unsigned n
)
{
UIntT2 t, b; // k, unsigned i, j; // counters
multiply_num_mn_unsigned( w, u, v, m, n );
// Теперь w[] содержит беззнаковое произведение.
// Корректируем вычитанием v*2**16m при и < 0 и
// вычитанием u*2**16n при v < Оif ((SIntT)u[m-1] < 0)
{
b = 0; // Инициализируем заемfor (j=0; j!=n; j++)
{
t = w[j+m] - v[j] - b;
w[j+m] = static_cast<UIntT>(t);
b = t >> (sizeof(UIntT2)*8 - 1); // 31;
}
}
if ((SIntT)v[n-1] < 0)
{
b = 0;
for(i=0; i!=m; i++)
{
t = w[i+n] - u[i] - b;
w[i+n] = static_cast<UIntT>(t);
b = t >> (sizeof(UIntT2)*8 - 1); // 31;
}
}
return;
}
inline void
multiply_num_mn_signed( UINT8 *w, const UINT8 *u, const UINT8 *v
, unsigned m, unsigned n )
{ multiply_num_mn_signed_impl<UINT8, INT8, UINT16>(w, u, v, m, n ); }
inline void
multiply_num_mn_signed( UINT16 *w, const UINT16 *u, const UINT16 *v
, unsigned m, unsigned n )
{ multiply_num_mn_signed_impl<UINT16, INT16, UINT32>(w, u, v, m, n ); }
inline void
multiply_num_mn_signed( UINT32 *w, const UINT32 *u, const UINT32 *v
, unsigned m, unsigned n )
{ multiply_num_mn_signed_impl<UINT32, INT32, UINT64>(w, u, v, m, n ); }
template<typename IntT, typename UIntT>
IntT multiply_impl(IntT u, IntT v, UIntT *pLowPart)
{
UIntT u0 = lo_part_s(u); // u & 0xFF;
IntT u1 = hi_part_s(u); // u >> 8;
UIntT v0 = lo_part_s(v); // v & 0xFF;
IntT v1 = hi_part_s(v); // v >> 8;
UIntT w0 = u0*v0;
IntT t = u1*v0 + hi_part_s(w0); // (w0 >> 8) ;
IntT w1 = lo_part_s(t); // t & 0xFF;
IntT w2 = hi_part_s(t); // t >> 8;
w1 = u0*v1 + w1;
if (pLowPart) *pLowPart = (UIntT)(w1 << calc_half_shift<UIntT>()) | (w0 & make_lo_half_mask<UIntT>() );
return u1*v1 + w2 + hi_part_s(w1); // (w1 >> 8) ;
}
template<typename IntT, typename UIntT>
IntT multiply_high_impl(IntT u, IntT v)
{
return multiply_impl(u, v, (UIntT*)0);
}
template<typename IntT>
IntT multiply_high( IntT a, IntT b)
{
return multiply_high_impl<UINT64, UINT64>( a, b );
}
template<> INT16 multiply_high<INT16> (INT16 a, INT16 b) { return multiply_high_impl<INT16, UINT16>( a, b ); }
template<> INT32 multiply_high<INT32> (INT32 a, INT32 b) { return multiply_high_impl<INT32, UINT32>( a, b ); }
template<> INT64 multiply_high<INT64> (INT64 a, INT64 b) { return multiply_high_impl<INT64, UINT64>( a, b ); }
template<> UINT16 multiply_high<UINT16>(UINT16 a, UINT16 b) { return multiply_high_impl<UINT16, UINT16>( a, b ); }
template<> UINT32 multiply_high<UINT32>(UINT32 a, UINT32 b) { return multiply_high_impl<UINT32, UINT32>( a, b ); }
template<> UINT64 multiply_high<UINT64>(UINT64 a, UINT64 b) { return multiply_high_impl<UINT64, UINT64>( a, b ); }
template <typename IntT>
IntT integer_pow( IntT i, unsigned p, bool *pOverflow = 0)
{
if (pOverflow) *pOverflow = false;
if (i==0 && p==0)
{
if (pOverflow) *pOverflow = true;
return 0;
}
if (!i) return 0;
if (!p) return 1;
IntT res = 1;
do{
IntT tmp = res*i;
IntT mh = multiply_high(res,i);
//sign_multiple( IntT a, IntT b )if (std::numeric_limits<IntT> :: is_signed )
{
if (i<0)
{
if (sign_multiple(tmp, res)>=0)
{
if (pOverflow) *pOverflow = true;
}
}
else
{
if (sign_multiple(tmp, res)<0)
{
if (pOverflow) *pOverflow = true;
}
}
int mhs = sign_multiple(res,i);
if ( (mhs<0 && mh!=(IntT)-1)
|| (mhs>0 && mh!=0)
)
{
if (pOverflow) *pOverflow = true;
}
}
else
{
if (mh>0)
{
if (pOverflow) *pOverflow = true;
}
}
res = tmp;
if (!res) { if (pOverflow) *pOverflow = true; return 0; }
--p;
} while(p);
return res;
}
template <typename IntT>
IntT integer_pow( IntT i, unsigned p, bool &bOverflow )
{
return integer_pow( i, p, &bOverflow );
}
Здравствуйте, Marty, Вы писали:
M> Понадобилось тут немного арифметики, решил сам наваять,
Понадобилось? Зачем? Почему не подошло готовое?
M> дело вообщем-то нехитрое ;)
Ну и раз нехитрое, то что же код деления не привёл тогда? )
M> Использовал эту
Ну, видимо невнимательно прочитал ту тему. Вот например:
M>// returns true if carry detected
M>inline bool signed_addition( INT32 a, INT32 b, INT32 *c = 0 )
M>{
M> INT32 r = a + b; if (c) *c = r; return !(b >= 0 == r >= a);
M>}
Переполнение для знаковых чисел — неопределённое поведение. Твоя проверка !(b >= 0 == r >= a) бессмысленна и неверна.
Компилятор, например, аналогично считает: http://goo.gl/2Gcbx4 — в результирующем кода функция всегда возвращает false. Ну и это очевидно почему так: если после сложения результат не вызвал переполнения, то возвращаем false. Если вызвал, то произошло UB и всё равно что функция будет возвращать, поэтому также возвращаем false, так как это не требует дополнительного кода. И оказывается, что от проверки ничего не зависит, и её можно смело игнорировать.
M>template<> bool integer_addition<INT16> (INT16 a, INT16 b, INT16 *c) { return signed_addition( a, b, c ); }
M>template<> bool integer_addition<INT16> (INT16 a, INT16 b, INT16 &c) { return signed_addition( a, b, c ); }
M>template<> bool integer_addition<INT32> (INT32 a, INT32 b, INT32 *c) { return signed_addition( a, b, c ); }
M>template<> bool integer_addition<INT32> (INT32 a, INT32 b, INT32 &c) { return signed_addition( a, b, c ); }
M>template<> bool integer_addition<INT64> (INT64 a, INT64 b, INT64 *c) { return signed_addition( a, b, c ); }
M>template<> bool integer_addition<INT64> (INT64 a, INT64 b, INT64 &c) { return signed_addition( a, b, c ); }
Кошмар. Шаблоны нужны чтобы избавиться от копипасты, а не чтобы сделать её в несколько раз больше. Ну вот зачем тебе тут специализация? Зачем вообще вызвать нешаблонные unsigned_addition из шаблонных integer_addition?
Я ещё понимаю, если бы ты разделял код для знаковых и беззнаковых типов (но даже для этого специализация не нужна), но вот это дублирование для всех размеров совсем ужасно.
Здравствуйте, watchmaker, Вы писали:
M>> Понадобилось тут немного арифметики, решил сам наваять,
W>Понадобилось? Зачем? Почему не подошло готовое?
Искать дольше
M>> дело вообщем-то нехитрое W>Ну и раз нехитрое, то что же код деления не привёл тогда? )
Типа подколол?
Там конечно несколько хитрее Пока не нужно, и возможно вообще не понадобится, не стал разбираться.
W>Ну, видимо невнимательно прочитал ту тему. Вот например:
Внимательно. Те компиляторы, которые я использую, все делают как нужно. Пока это устраивает.
W>
M>>template<> bool integer_addition<INT16> (INT16 a, INT16 b, INT16 *c) { return signed_addition( a, b, c ); }
W>
W>Кошмар. Шаблоны нужны чтобы избавиться от копипасты, а не чтобы сделать её в несколько раз больше. Ну вот зачем тебе тут специализация? Зачем вообще вызвать нешаблонные unsigned_addition из шаблонных integer_addition? W>Я ещё понимаю, если бы ты разделял код для знаковых и беззнаковых типов (но даже для этого специализация не нужна), но вот это дублирование для всех размеров совсем ужасно.
Так разделяю как раз. Хотя эта копипаста самому не очень нравится, потом подумаю, как переделать
W>
Здравствуйте, Marty, Вы писали:
W>>Понадобилось? Зачем? Почему не подошло готовое? M>Искать дольше
То есть вопрос правильности не стоит?
W>>Ну, видимо невнимательно прочитал ту тему. Вот например: M>Внимательно. Те компиляторы, которые я использую, все делают как нужно. Пока это устраивает.
Ясно. Хотя всё равно непонятно зачем оставлять явные ошибки, если можно сразу написать правильно.
M>Хотя эта копипаста самому не очень нравится,
В этом коде, кстати говоря, ещё один довольно сомнительный подход используется. Зачем нужно проводить вычисления над 8-,16-,32-х битными числами если длина машинного слова 64 бита (судя по наличию UINT64). Вместо этого просто сразу расширяешь все короткие числа до машинного слова, обрабатывать их как обычно, а потом, если нужно, при упаковке обратно проверять все переполнения.
Это же дико (и медленно) умножать два числа через полубайты, если у тебя регистр 64 или 32 бита. Ещё более дико становится если вспомнить, что много где есть возможность умножить пару машинных слов и получить результат двойной длины.
В общем, поэтому и спрашивал, зачем это нужно. Ибо сейчас код неправильный, небыстрый и некрасивый (сам просил попинать).
W>>Эта функция у тебя иногда возвращает отрицательный результат вследствие переполнения на минимальном целом. M>Это как? integer_abs должна возвращать неотрицательный результат, так? А вот integer_abs<int32_t>(-2147483648) возвращает число меньше нуля, хотя аргумент функции в int32_t помещается. Уж не знаю что ты с этим будешь делать, но явно тут прослеживается неконсистентность в подходе — одни функции пытаются сообщать об переполнениях, а другие их игнорируют.
Здравствуйте, watchmaker, Вы писали:
W>>>Понадобилось? Зачем? Почему не подошло готовое? M>>Искать дольше W>То есть вопрос правильности не стоит?
Стоит, конечно.
W>>>Ну, видимо невнимательно прочитал ту тему. Вот например: M>>Внимательно. Те компиляторы, которые я использую, все делают как нужно. Пока это устраивает. W>Ясно. Хотя всё равно непонятно зачем оставлять явные ошибки, если можно сразу написать правильно.
Тут даже не так. signed версии я добавил для ассортимента, копипастой с unsigned, и просто не подумал По идее, переполнения при сложении signed не может быть от слова никогда, а если оно вдруг как-то случилось, то это да, UB Со знаковыми может быть только займ при вычитании, надо будет подумать на эту тему.
M>>Хотя эта копипаста самому не очень нравится, W>В этом коде, кстати говоря, ещё один довольно сомнительный подход используется. Зачем нужно проводить вычисления над 8-,16-,32-х битными числами если длина машинного слова 64 бита (судя по наличию UINT64). Вместо этого просто сразу расширяешь все короткие числа до машинного слова, обрабатывать их как обычно, а потом, если нужно, при упаковке обратно проверять все переполнения.
Для ассортимента
Иногда пишу на TurboC 1.0 для 186 При нужде переделать под C можно без проблем, если все уже работает.
И на x86 маш. слово все-таки 32 разряда, возможно есть mul32, который дает 64х-битный результат, но не уверен, пусть с этим компилятор разбирается
Вопрос в том, как умножать числа, которые дают результат шире маш. слова. И для которых нет (или есть, но сильно нестандартные) интегральные типы с большей разрядностью.
W>Это же дико (и медленно) умножать два числа через полубайты, если у тебя регистр 64 или 32 бита. Ещё более дико становится если вспомнить, что много где есть возможность умножить пару машинных слов и получить результат двойной длины.
Полубайты?
Ты это уже к алгоритму умножения чисел произвольной длины ?
Где есть возможность умножения с расширением результата? На x86 вроде нет даже сейчас (64-128). Не знаю, есть ли на x64 умножение mul64 со 128-битным результатом. Вроде тоже нет.
Алгоритм умножения, который я стырил у Уоррена, как раз на это и расчитан.
W>В общем, поэтому и спрашивал, зачем это нужно. Ибо сейчас код неправильный, небыстрый и некрасивый (сам просил попинать).
За критику спасибо
W>>>Эта функция у тебя иногда возвращает отрицательный результат вследствие переполнения на минимальном целом. M>>Это как? W>integer_abs должна возвращать неотрицательный результат, так? А вот integer_abs<int32_t>(-2147483648) возвращает число меньше нуля, хотя аргумент функции в int32_t помещается.
А, понял мысль. Да, есть такая проблема. Когда вопрос встанет ребром, добавлю safe-версию.
W>Уж не знаю что ты с этим будешь делать, но явно тут прослеживается неконсистентность в подходе — одни функции пытаются сообщать об переполнениях, а другие их игнорируют.
Ну, те которые пытаются, они пытаются это явно сделать
Будет нужна safe integer_abs — сделаю отдельно, а в эту добавлю throw
Она вообще для ассортимента за компанию к integer_sign и sign_multiple сделана, даже и не используется пока
Я то что нужно было, сделал с большим внимание, а часть про запас накопипастил
W>>Это же дико (и медленно) умножать два числа через полубайты, если у тебя регистр 64 или 32 бита. Ещё более дико становится если вспомнить, что много где есть возможность умножить пару машинных слов и получить результат двойной длины.
M>Полубайты? M>Ты это уже к алгоритму умножения чисел произвольной длины ? M>Где есть возможность умножения с расширением результата? На x86 вроде нет даже сейчас (64-128). Не знаю, есть ли на x64 умножение mul64 со 128-битным результатом. Вроде тоже нет.
Вообще-то инструкция умножающая пару чисел и с результатом двойной длины была уже в самом первом 8086 процессоре. Соответственно, та же самая инструкция в x86-64 умножает пару 64-битных чисел с 128-битным результатом. То есть на 16-битном 8086 результатом будет 32-битное число, на 32-битном i386 результатом будет 64-битное число, и на современных x86-64 результатом умножения является 128-битное число.
Как раз наоборот, возможность перемножить пару чисел и получить результат такой же длины (то есть с отбрасыванием в случае переполнения) появилось в семействе x86 намного позже. А изначально все процессоры умножали с расширением (ну и продолжают это делать сейчас).
M>Алгоритм умножения, который я стырил у Уоррена, как раз на это и расчитан.
Ага. Поэтому если адаптировать его на x86, то нужно обрабатывать числа вдвое длиннее чем разрядность процессора — там как раз получается, что старшая/младшая половина помещается в регистр, а умножение пары регистров как раз с расширением и делается в процессоре.
W>>В этом коде, кстати говоря, ещё один довольно сомнительный подход используется. Зачем нужно проводить вычисления над 8-,16-,32-х битными числами если длина машинного слова 64 бита (судя по наличию UINT64). Вместо этого просто сразу расширяешь все короткие числа до машинного слова, обрабатывать их как обычно, а потом, если нужно, при упаковке обратно проверять все переполнения.
M>Для ассортимента ;) M>Иногда пишу на TurboC 1.0 для 186 ;) При нужде переделать под C можно без проблем, если все уже работает. M>И на x86 маш. слово все-таки 32 разряда, возможно есть mul32, который дает 64х-битный результат, но не уверен, пусть с этим компилятор разбирается ;) M>Вопрос в том, как умножать числа, которые дают результат шире маш. слова. И для которых нет (или есть, но сильно нестандартные) интегральные типы с большей разрядностью.
Ну я не много не про это. Ладно, пусть на процессоре нет операции умножения 32x32 -> 64, а есть только 32x32 -> 32 (с отбрасыванием старшей части) — не суть. И из-за этого ты используешь алгоритм Уоррена для 32-битных чисел. Но зачем ты его же используешь для 8-битных? Тут-то точно уже известно, что произведение двух 8-битных чисел влезет в 32 бита.
Ну то есть я понимаю, что он у тебя возник из-за того, что в шаблонах легко поменять параметр и получить реализацию для более коротких чисел. Но вот это-то как раз и плохо: вместо того чтобы умножать напрямую, у тебя используется менее эффективный алгоритм, созданный для совсем другой ситуации.
M>Я то что нужно было, сделал с большим внимание, а часть про запас накопипастил ;)
А тесты-то написал? Мне вот кажется, что код не должен работать. Например, вот эта строка t = u[i]*v[j] + w[i+j] + k;. Тут у тебя t, k имеют двойную длину, а аргументы у произведения — одинарную. При этом в языках C/C++ есть такая особенность — Integral promotion называется. Если, скажем, это соответственно 32- и 16-битные числа, то ошибка не проявляется, так как перед вычислением произведения происходит расширение аргументов до длины машинного слова, т.е. на твоём комптьютере скорее всего до 32-битного числа. А вот если это 64- и 32-битные числа, то автоматического расширения может и не произойти — и часть битов результата будет отброшена. Собственно, multiply_num_mn_unsigned_impl<UINT32, UINT64> у тебя в коде используется и непонятно как ты эту ошибку в тестах не поймал.
Здравствуйте, watchmaker, Вы писали:
W>Здравствуйте, Marty, Вы писали:
W>Вообще-то инструкция умножающая пару чисел и с результатом двойной длины была уже в самом первом 8086 процессоре. Соответственно, та же самая инструкция в x86-64 умножает пару 64-битных чисел с 128-битным результатом. То есть на 16-битном 8086 результатом будет 32-битное число, на 32-битном i386 результатом будет 64-битное число, и на современных x86-64 результатом умножения является 128-битное число.
Насчет 8086 уточнил — да, ты прав. Насчет 32 x86 — скорее всего тоже. Вот насчет 128х-разрядного результата не уверен, вроде совсем недавно где-то читал, что такого нет и только планируется.
А это собственно, и было нужно. Далее, если и есть — как это достать из C++? Писать на встроенном асме? Непереносимо. Писать две версии? Наверно, так. Или искать библиотеку, из которой можно быстренько дернуть кусочек кода? В пол-дня не уложился бы точно
В принципе, мне скорость тут абсолютно пофигу, я не криптографию ломать собрался Если приспичит, наверно подыщу замену.
M>>Алгоритм умножения, который я стырил у Уоррена, как раз на это и расчитан. W>Ага. Поэтому если адаптировать его на x86, то нужно обрабатывать числа вдвое длиннее чем разрядность процессора — там как раз получается, что старшая/младшая половина помещается в регистр, а умножение пары регистров как раз с расширением и делается в процессоре.
Да, я подзабыл, что у x86 есть, виновен Но опять же, как переносимо написать?
W>>>В этом коде, кстати говоря, ещё один довольно сомнительный подход используется. Зачем нужно проводить вычисления над 8-,16-,32-х битными числами если длина машинного слова 64 бита (судя по наличию UINT64). Вместо этого просто сразу расширяешь все короткие числа до машинного слова, обрабатывать их как обычно, а потом, если нужно, при упаковке обратно проверять все переполнения.
W>Ну я не много не про это. Ладно, пусть на процессоре нет операции умножения 32x32 -> 64, а есть только 32x32 -> 32 (с отбрасыванием старшей части) — не суть. И из-за этого ты используешь алгоритм Уоррена для 32-битных чисел. Но зачем ты его же используешь для 8-битных? Тут-то точно уже известно, что произведение двух 8-битных чисел влезет в 32 бита.
Я на этой версии тестировал С большими числами тесты сложнее сочинять.
W>Ну то есть я понимаю, что он у тебя возник из-за того, что в шаблонах легко поменять параметр и получить реализацию для более коротких чисел. Но вот это-то как раз и плохо: вместо того чтобы умножать напрямую, у тебя используется менее эффективный алгоритм, созданный для совсем другой ситуации.
Можно и убрать специализации, или переделать поэффективнее.
W>А тесты-то написал? Мне вот кажется, что код не должен работать. Например, вот эта строка t = u[i]*v[j] + w[i+j] + k;. Тут у тебя t, k имеют двойную длину, а аргументы у произведения — одинарную. При этом в языках C/C++ есть такая особенность — Integral promotion называется. Если, скажем, это соответственно 32- и 16-битные числа, то ошибка не проявляется, так как перед вычислением произведения происходит расширение аргументов до длины машинного слова, т.е. на твоём комптьютере скорее всего до 32-битного числа. А вот если это 64- и 32-битные числа, то автоматического расширения может и не произойти — и часть битов результата будет отброшена. Собственно, multiply_num_mn_unsigned_impl<UINT32, UINT64> у тебя в коде используется и непонятно как ты эту ошибку в тестах не поймал.
Надо посмотреть. Действительно, на максимуме разрядности мало потестил. Спасибо, что обратил внимание. Думаю, будет достаточно явно преобразовать к более широкому типу перед умножением.
M>Вот насчет 128х-разрядного результата не уверен, вроде совсем недавно где-то читал, что такого нет и только планируется.
Совсем недавно? Это есть в спецификации amd64, все процессоры, основанные на ней, уже более десятилетия умеют это. А сама спецификация даже ещё чуть старше.
M>Далее, если и есть — как это достать из C++? M>Но опять же, как переносимо написать?
Ну это как раз совсем не проблема:
#include <cstdint>
static inline void mul32(uint32_t a, uint32_t b, uint32_t& hi, uint32_t &lo) {
uint64_t t = a;
lo = t *= b;
hi = t >> 32;
}
Компиляторы (по крайней мере популярные gcc, clang и icc) знают, что uint64_t в данном коде — это всего лишь способ выразить желаемый результат, а не требование обязательно создать 64-битную переменную и как-то с ней работать даже на 32-битном процессоре (хотя они и это могут, для целей отладки, например). Собственно, если скомпилировать эту функцию на 32-битном процессоре x86, то она выродится в единственную инструкцию умножения двух 32-битных чисел, где 64-битный результат записывается в пару hi:lo.
Ну то есть код получается вполне переносим и при этом очень близок к тому что на самом деле происходит в процессоре.
Здравствуйте, watchmaker, Вы писали:
M>>Далее, если и есть — как это достать из C++? M>>Но опять же, как переносимо написать?
W>Ну это как раз совсем не проблема: W>...
Меня собственно, интересует такая версия:
static inline void mul128(uint64_t a, uint64_t b, uint64_t& hi, uint64_t &lo) {
uint128_t t = a;
lo = t *= b;
hi = t >> 64;
}
Конечно же, для умножения 32x32 я бы не стал огород городить
Здравствуйте, watchmaker, Вы писали:
W>А тесты-то написал? Мне вот кажется, что код не должен работать. Например, вот эта строка t = u[i]*v[j] + w[i+j] + k;. W>Тут у тебя t, k имеют двойную длину, а аргументы у произведения — одинарную. При этом в языках C/C++ есть такая особенность — Integral promotion называется.
Спасибо за замечание. Действительно был баг в моей версии. static_cast до широкого числа вылечил.
64x64 -> 128 — теперь работает как нужно
Здравствуйте, Marty, Вы писали:
M>Меня собственно, интересует такая версия: M>
static inline void mul128(uint64_t a, uint64_t b, uint64_t& hi, uint64_t &lo) {
M> uint128_t t = a;
M> lo = t *= b;
M> hi = t >> 64;
M>}
Хорошо, так что тут конкретно интересует?
Под архитектуру x86-64 это пример компилируется gcc, clang, icc. После компиляции — та же единственная инструкция умножения, только имена регистров другие.
Или ты намекаешь, что раз в текущем стандарте языка не упомянут тип uint256_t, то такой подход не масштабируется на гипотетические 128-битные процессоры? Ну так так всегда будет — сначала вводится платформо-зависимое расширение (в данном случае типы вроде __uint256_t или __uint128_t), а потом уже это попадает в стандарт (да и то, только если достаточно широко начинает использоваться).
Всё-таки языки C/C++ тем и отличаются, что в стандарте не описывается единственно верное поведение, в вместо этого существует, например, implementation defined behaviour. Это, конечно, не добавляет переносимости программам, но уж такой язык. В С++ ведь не только длины целых чисел меняются, так даже two's complement представление знаковых не гарантируется. Да даже размер байта не всегда 8 бит (и так действительно бывает). С этой точки зрения у тебя в программу уже содержится явно непереносимые операции вроде mask <<= 8;. Хочешь определённости и переносимости — пиши на какой-нибудь Java, а не на C++. А на С/C++ тебе самому придётся задумываться о том, является ли, например, long int 32-битным или 64-битным числом, и что твоя программа будет делать если таки не является.
Здравствуйте, watchmaker, Вы писали:
M>>Меня собственно, интересует такая версия: M>>
static inline void mul128(uint64_t a, uint64_t b, uint64_t& hi, uint64_t &lo) {
M>> uint128_t t = a;
M>> lo = t *= b;
M>> hi = t >> 64;
M>>}
W>Хорошо, так что тут конкретно интересует? W>Под архитектуру x86-64 это пример компилируется gcc, clang, icc. После компиляции — та же единственная инструкция умножения, только имена регистров другие.
Мне бы для MSVC 2005
W>Или ты намекаешь, что раз в текущем стандарте языка не упомянут тип uint256_t, то такой подход не масштабируется на гипотетические 128-битные процессоры?
А что, в текущем стандарте uint128_t уже есть?
W>Ну так так всегда будет — сначала вводится платформо-зависимое расширение (в данном случае типы вроде __uint256_t или __uint128_t), а потом уже это попадает в стандарт (да и то, только если достаточно широко начинает использоваться).
Было бы хорошо в стандарте иметь функции/операции для арифметики, с проверкой переполнения. Во всех железных архитектурах так или иначе какие-то флаги выставляются, а из C++, да и из C, до них не добраться. Большой недостаток для языка, который декларирует поддержку самых низкоуровневых деталей.Имхо, конечно. "C", вон, вообще как "переносимый асм" декларируется, а таких нужных фич нет.
W>Всё-таки языки C/C++ тем и отличаются, что в стандарте не описывается единственно верное поведение, в вместо этого существует, например, implementation defined behaviour. Это, конечно, не добавляет переносимости программам, но уж такой язык. В С++ ведь не только длины целых чисел меняются, так даже two's complement представление знаковых не гарантируется. Да даже размер байта не всегда 8 бит (и так действительно бывает). С этой точки зрения у тебя в программу уже содержится явно непереносимые операции вроде mask <<= 8;. Хочешь определённости и переносимости — пиши на какой-нибудь Java, а не на C++. А на С/C++ тебе самому придётся задумываться о том, является ли, например, long int 32-битным или 64-битным числом, и что твоя программа будет делать если таки не является.
О первом выделенном — помню, лень было искать, как называется количество бит на байт в limits
По второму — мне особо не интересно, какие там у чего размеры реально. Мне интересно отслеживать факты переполнения при операциях с моими данными.
А что, в Java такие вещи можно отслеживать? Или там исключениями будут кидаться?
Здравствуйте, Marty, Вы писали:
W>>Или ты намекаешь, что раз в текущем стандарте языка не упомянут тип uint256_t, то такой подход не масштабируется на гипотетические 128-битные процессоры? M>А что, в текущем стандарте uint128_t уже есть?
В текущем стандарте даже uint32_t помечен как optional . Это к тому, что с целыми числами в стандарте языка всё совсем не так просто. Тут именно, либо ты используешь что-то вроде int или size_t, имеющие некие гарантии на диапазон, либо вступаешь в область платформо-зависимого поведения. Ну, понятно, что в большинстве случаев int и size_t весьма неплохой выбор и особо можно не переживать. И, отвечая на твой вопрос, uint128_t в текущем стандарте нет. Но, опять же, это довольно мало значит.
M> Во всех железных архитектурах так или иначе какие-то флаги выставляются, а из C++, да и из C, до них не добраться. Большой недостаток для языка, который декларирует поддержку самых низкоуровневых деталей.Имхо, конечно. "C", вон, вообще как "переносимый асм" декларируется, а таких нужных фич нет.
«Переносимый асм» не означает, что ты скопировал программу на другую архитектуру и она заработала. Это означает лишь что после копирования тебе придётся поправить не всю программу (как это было бы с ассемблером), а лишь небольшую её часть. При этом язык сразу предоставляет готовые примитивы вроде sizeof(void*), которые очевидно меняются часто. Но он никак не поможет изменить взаимодействие файловой системой, с ОС, права доступа к памяти, особенности ввода-вывода — тут всё настолько разнообразно, что просто нельзя это никак уровнять.
Часто стараются засунуть это в какие-то библиотеки с глаз долой. Но хотя библиотека, реализующая putchar('a'), и написана на переносимом С, но вот именно печать символа для каждой платформы написана по своему непереносимо.
Аналогично и с сигнализацией переполнения — везде всё настолько по-разному реализовано, что какого-то единого интерфейса предложить нельзя. В лучшем случае что может тут предложить С/C++ — это сноску «implementation defined behaviour» (а часто всё даже заканчивается «undefined behaviour»). Хочешь сигнализации (и прочего детерминизма) — смотри на расширения компилятора — обычно там уже много полезных примитивов есть.
M>А что, в Java такие вещи можно отслеживать? Или там исключениями будут кидаться?
А там виртуальная машина жёстко зафиксирована. Неважно какой у тебя процессор, архитектура, интерпретатор или jit-компилятор, — изволь эмулировать то что написано. Например, тот же long в java с точки зрения программы всегда 64-битный и в two's complement представлении — это гарантируется. Если процессор так не умеет, то реализация обязана, например, программной эмуляцией добиться соответствующего поведения.
Здравствуйте, watchmaker, Вы писали:
W>Здравствуйте, Marty, Вы писали:
W>>>Или ты намекаешь, что раз в текущем стандарте языка не упомянут тип uint256_t, то такой подход не масштабируется на гипотетические 128-битные процессоры? M>>А что, в текущем стандарте uint128_t уже есть? W>В текущем стандарте даже uint32_t помечен как optional .
Если на платформе нет int32 — придется обходится без него ;( Но потребность в арифметике, которая превышает разрядность аппаратуры, никуда не исчезает.
W>Это к тому, что с целыми числами в стандарте языка всё совсем не так просто. Тут именно, либо ты используешь что-то вроде int или size_t, имеющие некие гарантии на диапазон, либо вступаешь в область платформо-зависимого поведения. Ну, понятно, что в большинстве случаев int и size_t весьма неплохой выбор и особо можно не переживать. И, отвечая на твой вопрос, uint128_t в текущем стандарте нет. Но, опять же, это довольно мало значит.
Ну, sizeof + CHAR_BITS помогают
M>> Во всех железных архитектурах так или иначе какие-то флаги выставляются, а из C++, да и из C, до них не добраться. Большой недостаток для языка, который декларирует поддержку самых низкоуровневых деталей.Имхо, конечно. "C", вон, вообще как "переносимый асм" декларируется, а таких нужных фич нет. W>«Переносимый асм» не означает, что ты скопировал программу на другую архитектуру и она заработала.
Кто же спорит
W>Но он никак не поможет изменить взаимодействие файловой системой, с ОС, права доступа к памяти, особенности ввода-вывода — тут всё настолько разнообразно, что просто нельзя это никак уровнять.
Ну, к аппаратной платформе это мало отношения имеет.
W>Часто стараются засунуть это в какие-то библиотеки с глаз долой. Но хотя библиотека, реализующая putchar('a'), и написана на переносимом С, но вот именно печать символа для каждой платформы написана по своему непереносимо.
Ну и хотелось бы подобных стандартных функций, которые реализованы максимально эффективно для конкретной платформы, и поддерживаются производителем, а не группой опен-сорц энтузиастов или вообще своими силами.
W>Аналогично и с сигнализацией переполнения — везде всё настолько по-разному реализовано, что какого-то единого интерфейса предложить нельзя.
Есть процессоры без флаговых регистров? А язык C для них есть?
Intel, ARM, MIPS — есть флаги; контроллеры всякие — все, что видел, тоже с флаговым регистром в том или ином виде. Поведение несколько отличается, но это как раз повод просто упрятать конкретную реализацию в стд библиотеку. Просто стандарт обязывал бы производителя платформы поддерживать еще пару-тройку функций, вот и все.
W>В лучшем случае что может тут предложить С/C++ — это сноску «implementation defined behaviour» (а часто всё даже заканчивается «undefined behaviour»). Хочешь сигнализации (и прочего детерминизма) — смотри на расширения компилятора — обычно там уже много полезных примитивов есть.
M>>А что, в Java такие вещи можно отслеживать? Или там исключениями будут кидаться? W>А там виртуальная машина жёстко зафиксирована. Неважно какой у тебя процессор, архитектура, интерпретатор или jit-компилятор, — изволь эмулировать то что написано. Например, тот же long в java с точки зрения программы всегда 64-битный и в two's complement представлении — это гарантируется. Если процессор так не умеет, то реализация обязана, например, программной эмуляцией добиться соответствующего поведения.
Ну, это да. Но как там с переполнениями?
Ладно, мы тут уже далеко отошли от изначальной темы, надо завязывать треп
Здравствуйте, Кодт, Вы писали:
M>> Понадобилось тут немного арифметики, решил сам наваять, дело вообщем-то нехитрое
К>Не смотрел ли на GMP? http://en.wikipedia.org/wiki/GNU_Multi-Precision_Library К>или на Boost::Multi-Precision?
По поводу GMP — о ней слышал. Но 1) GNU, 2) Много за собой тянет. Мне надо просто пару функций в хидер-онли варианте.
По поводу Boost::Multi-Precision — буст в своих проектах стараюсь не использовать. Если есть указание свыше, тогда — да.
Причины просты: буст затягивает, чем больше его используешь, тем больше его использовать приходится. Я пробовал им пользоваться, но при обычно, если что-то нужно, оказывается, что в нынешней версии буста нужная либа содержит баги или функционал недоработан до нужного тебе, переезжаешь на новую версию, и тут начинают вылезать всяческие проблемы, связанные с тем, что тут устарело, там интерфейс поменяли, это вообще выкинули.
Я просто жду, когда основные нужности переедут в стандарт C++ и будут поддерживаться всеми основными компиляторами в одинаково полной мере
Вообще, конечно, если по арифметике будут еще задачи всплывать, то буду подыскивать что-то на замену велосипеду.
Нашел тут еще чей-то "мопед" — cpp-bigint. Пока не смотрел детально, но вроде довольно легковесный — по одному cpp и h файлу. Из буста использует только boost/static_assert.hpp, да и то если указать явно.
Здравствуйте, Marty, Вы писали:
W>>Аналогично и с сигнализацией переполнения — везде всё настолько по-разному реализовано, что какого-то единого интерфейса предложить нельзя. M>Есть процессоры без флаговых регистров? А язык C для них есть? M>Intel, ARM, MIPS — есть флаги;
Ты явно что-то путаешь. У MIPS флагового регистра в привычном виде нет. У него есть команды перехода по условию и установки целевого регистра по условию (в 0 или 1), а условие выглядит как сравнение src1 и src2, при этом src1 — регистр, а src2 — регистр или константа. Условия — стандартный набор типа "равно/больше/меньше" с вариантами со знаком или без. Такая же ситуация на Alpha. При разработке обоих тут следовали идее, что флаговый регистр (PSW) — узкое место и не надо завязываться на него для всех операций.
Для примера возьмём MIPS код в libgmp (mpn/mips32/add_n.asm, сложение длинных натуральных), комментирую:
C INPUT PARAMETERS
C res_ptr $4
C s1_ptr $5
C s2_ptr $6
C size $7
ASM_START()
PROLOGUE(mpn_add_n)
lw $10,0($5); младшее слово первого аргумента
lw $11,0($6); младшее слово второго аргумента
addiu $7,$7,-1
and $9,$7,4-1 C number of limbs in first loop
beq $9,$0,.L0 C if multiple of 4 limbs, skip first loop
move $2,$0; инициализация переноса в 0
subu $7,$7,$9
; входим со значениями в текущей позиции, уже загруженными в $10 и $11
.Loop0: addiu $9,$9,-1; декремент переменной цикла - позиции в массиве
lw $12,4($5); прочли следующее слово из 1-го аргумента
addu $11,$11,$2; прибавили перенос к текущему значению из 2-го аргумента
lw $13,4($6); прочли следующее слово из 2-го аргумента
sltu $8,$11,$2; запись результата сравнения на меньше без знака,
; то есть $8 получит 1 только если перенос с предыдущего цикла
; дал беззнаковое переполнение, иначе ставится в 0
addu $11,$10,$11; складываем два входных значения текущей позиции
sltu $2,$11,$10; запись результата сравнения на меньше без знака, аналогично,
; запись факта переполнения
sw $11,0($4); выгрузка суммы на текущем шаге в массив результата
or $2,$2,$8; объединение двух признаков переноса в один логическим "или"
addiu $5,$5,4; инкремент указателя 1-го аргумента на 1 слово
addiu $6,$6,4; инкремент указателя 2-го аргумента на 1 слово
move $10,$12; следующее прочтённое слово из 1-го аргумента сделать текущим
move $11,$13; следующее прочтённое слово из 2-го аргумента сделать текущим
bne $9,$0,.Loop0
addiu $4,$4,4
Дальнейший код (ещё один цикл, более крупными прыжками) пропустил. Я надеюсь, по комментариям видно, как получается сумма и вычисление переноса. На C один шаг такого суммирования внутри цикла можно было бы записать так:
Громоздко? Да. Дольше варианта, когда явно ставится и анализируется флаг переноса?
Да, если считать на количество команд.
Но: работает? Да, работает, корректно. И в описанном сишном варианте гарантированно работает независимо от архитектуры и без ассемблерных вставок (которые, конечно, на каком-нибудь x86 ускорят итерацию в несколько раз).
Аналогичные проблемы, но несколько более облегчённые, были на S/360 (до введения в S/390 команд типа ALC), несмотря на то, что это CISC. Там надо было точно так же сделать два сложения, от обоих сохранить выходной бит 1 поля CC и сделать логический OR для получения признака переноса на следующий цикл. По сравнению с этим привычный тебе формат (появившийся впервые, насколько я смог раскопать, в PDP-11), это уже позднее сверхудачное нововведение.
M> контроллеры всякие — все, что видел, тоже с флаговым регистром в том или ином виде. Поведение несколько отличается, но это как раз повод просто упрятать конкретную реализацию в стд библиотеку. Просто стандарт обязывал бы производителя платформы поддерживать еще пару-тройку функций, вот и все.
Как видишь, нет, не с флаговым регистром. С другой стороны, я полностью согласен с тем, что функцию, идентичную по поведению такой sum_step, имеет смысл подумать внести в стандарт, хотя бы в необязательные, но в расширения со стандартизованным интерфейсом. Ну а на каком-нибудь x86 она превратится в банальную цепочку из shr + adc + setc, с возможностью устранения крайних элементов в ней, если между вызовами CF не меняется.
Операции с числами со знаком тоже легко делаются через такое же, просто старший бит старшего слова значения определяет знак Но обрати внимание, что тот же GMP сделал mpn нижним уровнем, а знаки чисел уже над ним. При "реально" длинной арифметике так эффективнее.
M>>>А что, в Java такие вещи можно отслеживать? Или там исключениями будут кидаться? W>>А там виртуальная машина жёстко зафиксирована. Неважно какой у тебя процессор, архитектура, интерпретатор или jit-компилятор, — изволь эмулировать то что написано. Например, тот же long в java с точки зрения программы всегда 64-битный и в two's complement представлении — это гарантируется. Если процессор так не умеет, то реализация обязана, например, программной эмуляцией добиться соответствующего поведения. M>Ну, это да. Но как там с переполнениями?
Не проверяются, там арифметика гарантированно кольцевая по модулю (2**32 или 2**64).
Зато ты можешь проверить результаты как минимум сложения или вычитания практически таким же способом, как я показал: если c = a+b, то проверка на переполнение делается через:
((b >= 0) ? (c < a) : (c > a))
Для GCC аналогичное можно получить, скомпилировав с -fwrapv (для целого входного файла для версий до 4.4, с грануляцией до функции, начиная с 4.4). Clang знает тот же -fwrapv для входного файла.
Вот с проверкой умножения, конечно, сложнее всего. Обратное деление слишком дорого. Если нет ассемблерной поддержки, лучше брать входные значения размером в половину доступного слова. Будет дольше, зато надёжно.
Здравствуйте, Marty, Вы писали:
M>inline bool unsigned_addition( UINT16 a, UINT16 b, UINT16 *c = 0 ) M>{ M> UINT16 r = a + b; if (c) *c = r; return !(r>=a && r>=b); M>}
Я правильно понимаю, что возвращается во всех функциях 1, если нет ошибок, и 0, если переполнение?
Если так, то часть проверок тут излишня — для детекта переполнения достаточно одного из условий r < a или r < b. В случае переполнения они оба будут нарушены одновременно.
M>inline bool signed_addition( INT64 a, INT64 b, INT64 *c = 0 ) M>{ M> INT64 r = a + b; if (c) *c = r; return !(b >= 0 == r >= a); M>}
, компиляторы нынче пошли очень подлые и этот стиль проверки чреват боком (™). Написав сложение без условия, ты подсказываешь, что переполнения не будет. Варианты лечения:
1. Заставить скомпилировать этот кусок с -fno-strict-overflow или даже -fwrapv (gcc и последователи).
2. Поставить проверку вначале:
inline bool signed_addition( INT64 a, INT64 b, INT64 *c = 0 )
{
if (b > 0 && a > INT64_MAX - b)
return 0;
if (b < 0 && a + b < INT64_MIN)
return 0;
if (c)
*c = a + b;
return 1;
}
Тут у меня чуть другая семантика (при переполнении нет присвоения). Если нужно идеальное совпадение с предыдущей — придётся ещё больше извращаться, см. вариант 3.
3. Пересчитать в unsigned и сравнить флаги самому:
inline bool signed_addition( INT64 a, INT64 b, INT64 *c = 0 )
{
UINT64 mask = 0x8000000000000000ull;
UINT64 r = (UINT64) a + (UINT64)b;
if (c) *c = (INT64)r;
return (a&mask) != (b&mask) || (r&mask) == (a&mask);
}
Аналогично для вычитания, с симметричным переворотом подходов.
M>template<typename IntT> M>int sign_multiple( IntT a, IntT b ) M> { return integer_sign(a)==integer_sign(b) ? 1 : -1; }
Умножить на CHAR_BIT из <limits.h> и результат умножения распополамить. (Если нечётное — вылететь с громким криком.)
M>template<typename IntT> M>IntT hi_part_s( IntT i ) M> { return i >> calc_half_shift<IntT>(); }
Вот тут есть грабелька: для знакового значения не гарантируется конкретное исполнение сдвига — с размножением знака или нет. Поэтому hi_part_s<int32_t>(-65536) может оказаться и -1, и 65535.
M>template<typename UIntT, typename UIntT2> M>void multiply_num_mn_unsigned_impl( UIntT *w M> , const UIntT *u, const UIntT *v M> , unsigned m, unsigned n M> ) M>{ M> UIntT2 k, t; //, b; M> unsigned i, j; // counters
M> for (i=0; i<m; i++) M> w[i] = 0;
M> for(j=0; j!=n; j++) M> { M> k = 0; M> for(i=0; i!=m; i++) M> { M> t = u[i]*v[j] + w[i+j] + k;
Если на обычных платформах (где int == int32_t), то для пары <INT16, INT32> это будет работать, а вот для пары <INT32, INT64> уже нет — ты забыл проконвертировать u[i] и v[j] в более высокую разрядность (INT64 в последнем случае) до их умножения.
Умножение двух 32-разрядных чисел не будет расширяться до 64 даже в случае переполнения, это справедливо даже на x86-64.
Если на чём-то более мелкое с
M> w[i+j] = static_cast<UIntT>(lo_part_s(t)); // t; // (Т.е., t & OxFFFF).
Мнэээ. Во-первых, нафига тебе lo_part_s(), когда lo_part() и делает то же самое, и яснее по названию? Во-вторых, а почему lo_part() вообще работает со знаковыми, а не беззнаковыми?
M> k = hi_part_s(t); // t >> 16;
Аналогично — тут лучше hi_part(). Кстати, если вернуться к их определению:
M>template<typename IntT> M>IntT lo_part( IntT i ) { return i & make_lo_half_mask<IntT>(); }
M>template<typename IntT> M>IntT hi_part( IntT i ) M> { return (i >> calc_half_shift<IntT>()) & make_lo_half_mask<IntT>(); }
если бы ты сразу написал как беззнаковые, финальный and с маской был бы не нужен. (Опять-таки, считая, что количество бит в слове чётное. Надеюсь, на M-20 с её 45-битным словом ты это переносить не будешь.)
M>template<typename UIntT, typename SIntT, typename UIntT2> M>void multiply_num_mn_signed_impl( UIntT *w M> , const UIntT *u, const UIntT *v M> , unsigned m, unsigned n M> ) M>{ M> UIntT2 t, b; // k, M> unsigned i, j; // counters M> multiply_num_mn_unsigned( w, u, v, m, n );
M> // Теперь w[] содержит беззнаковое произведение. M> // Корректируем вычитанием v*2**16m при и < 0 и
Заметно, какие у тебя шрифты на экране Или это из Уоррена?
s/и/u/
M> // вычитанием u*2**16n при v < О M> if ((SIntT)u[m-1] < 0)
Кстати, почему у тебя местами static_cast, а местами конверсия в стиле C?
M> { M> b = 0; // Инициализируем заем M> for (j=0; j!=n; j++) M> { M> t = w[j+m] — v[j] — b; M> w[j+m] = static_cast<UIntT>(t); M> b = t >> (sizeof(UIntT2)*8 — 1); // 31; M> } M> }
А почему не используешь свои же unsigned_subtraction?
M>template<typename IntT, typename UIntT> M>IntT multiply_high_impl(IntT u, IntT v) M>{ M> return multiply_impl(u, v, (UIntT*)0); M>}
(сам просил попинать).
2005 
Перекуём баги на фичи!
