Qraizer
Advanced Member | Редактировать | Профиль | Сообщение | ICQ | Цитировать | Сообщить модератору
/***************************************************************************\
|************** Решение квадратного уравнения в **************|
|******************** во время компиляции ********************|
\***************************************************************************/
#include <iostream>
#include <limits>
#include <cmath>
#include <cstddef>
#ifdef SHOW_DEBUG_INFO
/*
Отладочный шаблон для вывода значения I типа T в виде сообщения компилятора
об ошибке, в котором присутствует значение I.
Вложенный шаблон используется для более лёгкого нахождения требуемых данных
в тексте сообщения. Так значение I должно быть выведено в самом начале,
в противном случае оно будет закопано в глубине многострочного сообщения.
По крайней так происходит с Intel C++ Compiler-ом. Впрочем, в связи с тем,
стандарт не определяет - даже формально - форматы текстов об ошибках, нет
никакой гарантии, что информация будет читабельной. Например, Visual C++ 7.1
вообще выводит значение I в hex.
Значение выводится в compile-time, да ещё и в виде текста об ошибке, поэтому
при этом компиляция закончится неудачей.
*/
template <typename T, T I>
class outT
{
template <typename valueT, valueT J> struct OUTPUT;
public:
OUTPUT<T, I> _dummy;
};
#endif
// Шаблон класса - делителя двух чисел с математическим округлением
template <typename T, T X, T Y>
class Div
{
static const T quot=X/Y, // частное
part=quot-(quot<0 && X!=quot*Y),// поправка для отрицательных
Z1 =part*Y, // точность с недостатком
Z2 =(part+1)*Y; // точность с избытком
public:
static const T value=X-Z1 > Z2-X ? part+1 : // выбрать лучшее
X-Z1 !=Z2-X ? part :
part%2 == 0 ? part : part+1; // чётное
};
/*
К сожалению, форматы с плавающей точкой в качестве аргументов
шаблонов не разрешаются. Поэтому будет использоваться формат
с фиксированной точкой. Даже long на 32-битных платформах черезчур
мал для достаточной точности, поэтому тип аргументов параметризируется.
По это причине в классах вместо enum используются static const.
*/
// Параметры: тип аргументов, коэффициенты и количество разрядов
// в дробной части (по умолчанию - точность float)
template <
typename T,
T Ax2, T Bx, T C,
int Frac=6
>
class Qwur
{
// Вычисление корней по коэффициентам при неизвестном, дискриминанте и
// поправочном множителе вида std::pow(10, Frac)
template <
typename valueT,
valueT A, valueT B, valueT Discr,
valueT Pow
>
class Roots
{
/*
По хорошему нужна ещё специализация для нуля, но так как здесь она
не используется, то я её и не определил.
*/
// Вычисление квадратного корня методом Ньютона.
// Параметры: подкоренное значение и поправочный множитель.
template <
typename valueT1,
valueT1 X0,
valueT1 Pow1
>
class Sqrt
{
// Итерации по Ньютону (объявление).
// Параметры: подкоренное значение, очередное приближение
// и поправочный множитель.
template <
typename valueT2,
valueT2 X, valueT2 Xn,
valueT2 Pow2
>
struct Iterate;
// Выбор между продолжением итераций и их завершением.
// Параметры: подкоренное значение, предыдущее и очередное приближения,
// поправочный множитель и признак перехода итерации через точное значение.
// Последний параметр важен для частичной специализации.
template <
typename valueT2,
valueT2 X, valueT2 Xpred, valueT2 XN,
valueT2 Pow2,
bool Gt
>
struct Compare
{
static const valueT2 value=Iterate<valueT2, X, XN, Pow2>::value;
};
// Специализация для окончания рекурсии.
template <
typename valueT2,
valueT2 X, valueT2 XN,
valueT2 Pow2
>
struct Compare<valueT2, X, XN, XN, Pow2, false>
{ // Округление аналогично Div<>
static const valueT2 value=X-XN*XN > (XN+1)*(XN+1)-X ? XN+1 :
X-XN*XN !=(XN+1)*(XN+1)-X ? XN :
XN%2 ==0 ? XN : XN+1;
};
// Специализация для предотвращения бесконечной рекурсии в случае, если две
// (а может и больше) итерации начинают чередоваться.
// В целочисленной арифметике это нередкий случай.
template <
typename valueT2,
valueT2 X, valueT2 Xpred, valueT2 XN,
valueT2 Pow2
>
struct Compare<valueT2, X, Xpred, XN, Pow2, false>
{ // Округление аналогично Div<>, но...
static const valueT2 value=X-Xpred*Xpred > XN*XN-X ? XN :
X-Xpred*Xpred !=XN*XN-X ? Xpred :
// ... учитывая, что Xpred и XN могут отличаться более чем на 1,
// берём их среднее. Если всё же ровно на 1, то оно уже будет
// округлено до чётного, иначе - не факт.
Div<valueT2, XN+Xpred, 2>::value%2 == 0 ?
Div<valueT2, XN+Xpred, 2>::value : Div<valueT2, XN+Xpred, 2>::value+1;
};
// Итерации по Ньютону (определение).
template <
typename valueT2,
valueT2 X, valueT2 Xn,
valueT2 Pow2
>
struct Iterate // X = /X +A/X \ / новое
{ // n+1 \ n / n/ / 2 приближение
static const valueT2 Xnext=Div<valueT2, Xn+Div<valueT2, X*Pow2, Xn>::value, 2>::value;
public: // Определиться с окончанием/продолжением итераций
static const valueT2 value=Compare<valueT2, X, Xn, Xnext, Pow2, (Xn > Xnext)>::value;
};
public:
// Запустить итерации. Начальное приближение - подкоренное значение.
static const valueT1 value=Iterate<valueT1, X0, X0, Pow1>::value;
};
// Шаблон вычисления sign(x)
template <typename valueT1, valueT1 X>
struct Sign
{
static const valueT1 value = X<0 ? -1 : X==0 ? 0 : +1;
};
// Объявление шаблона выбора алгоритма вычисления корней
// в зависимости от знака дискриминанта.
// Параметры: коэффициенты при неизвестном, дискриминант,
// поправочный множитель и знак дискриминанта.
template <
typename valueT1,
valueT1 A, valueT1 B, valueT1 D,
valueT1 P,
int
>
struct Select;
// Специализация для отрицательного дискриминанта
template <
typename valueT1,
valueT1 A, valueT1 B, valueT1 D,
valueT1 P
>
class Select<valueT1, A, B, D, P, -1>
{
static const valueT1 sqrtD=Sqrt<valueT1, -D, P>::value;
public:
static const valueT1 x1Re =Div<valueT1, -B*P , 2*A>::value,
x2Re =x1Re ,
x2Im =Div<valueT1, sqrtD*P, 2*A>::value,
x1Im =-x2Im;
#ifdef SHOW_DEBUG_INFO
outT<valueT, D> _dummy1;
outT<valueT, sqrtD> _dummy2;
#endif
};
// Специализация для нулевого дискриминанта (не обязательна, но с ней проще)
template <
typename valueT1,
valueT1 A, valueT1 B, valueT1 D,
valueT1 P
>
struct Select<valueT1, A, B, D, P, 0>
{
static const valueT1 x1Re =Div<valueT1, -B*P, 2*A>::value,
x2Re =x1Re ,
x1Im =0 ,
x2Im =0;
#ifdef SHOW_DEBUG_INFO
outT<valueT, D> _dummy1;
#endif
};
// Специализация для положительного дискриминанта
template <
typename valueT1,
valueT1 A, valueT1 B, valueT1 D,
valueT1 P
>
class Select<valueT1, A, B, D, P, +1>
{
static const valueT1 sqrtD=Sqrt<valueT1, D, P>::value;
public:
static const valueT1 x1Re =Div<valueT1, (-B-sqrtD)*P, 2*A>::value,
x2Re =Div<valueT1, (-B+sqrtD)*P, 2*A>::value,
x1Im =0 ,
x2Im =0;
#ifdef SHOW_DEBUG_INFO
outT<valueT, D> _dummy1;
outT<valueT, sqrtD> _dummy2;
#endif
};
// Получить результаты и вытащить наружу
public:
static const valueT
x1Re=Select<valueT, A, B, Discr, Pow, Sign<valueT, Discr>::value>::x1Re,
x2Re=Select<valueT, A, B, Discr, Pow, Sign<valueT, Discr>::value>::x2Re,
x1Im=Select<valueT, A, B, Discr, Pow, Sign<valueT, Discr>::value>::x1Im,
x2Im=Select<valueT, A, B, Discr, Pow, Sign<valueT, Discr>::value>::x2Im;
};
// Вычисление целой степени десяти для преобразования количества дробных
// разрядов в значении с фиксированной точкой в поправочный множитель
template <T I> struct Pow10 { static const T value = Pow10<I-1>::value * 10; };
template <> struct Pow10<0>{ static const T value = 1; };
// Вычисление дискриминанта
template <
typename valueT,
valueT A, valueT B, valueT C,
valueT Pow
>
struct Discr
{
static const valueT value=Div<valueT, B*B-4*A*C, Pow>::value;
};
static const T pow =Pow10<Frac>::value; // поправочный множитель
static const T discr=Discr<T, Ax2, Bx, C, pow>::value;// дискриминант
public:
// Найти корни и вытащить наружу - x1=a1+i*b1, x2=a2+i*b2
static const T x1Re=Roots<T, Ax2, Bx, discr, pow>::x1Re, // a1
x2Re=Roots<T, Ax2, Bx, discr, pow>::x2Re, // a2
x1Im=Roots<T, Ax2, Bx, discr, pow>::x1Im, // b1
x2Im=Roots<T, Ax2, Bx, discr, pow>::x2Im; // b2
};
/*
Для тестирования правильности работы шаблона Qwur<> проверяем результаты
вычислениями в run-time и с плавающей точкой
*/
// Простой манипулятор для вывода табуляции
template <typename Ch, class Tr>
std::basic_ostream<Ch, Tr>& tab(std::basic_ostream<Ch, Tr>& os)
{
// static const Ch tabSymbol=os.widen('\t');
// return os.width(10)<<tabSymbol;
os.width(10);
return os;
}
// Получить коэффициенты и вернуть составные части комплексных корней
void qwur(float a, float b, float c, float& x1re, float& x1im, float& x2re, float& x2im)
{
float d=b*b-4*a*c;
// Отладочный вывод дискриминанта и квадратного корня его модуля
std::cout<<"d="<<d<<'\t'<<"sqrt(|d|)="<<std::sqrt(abs(d))<<std::endl;
// Разумная погрешность при сравнении на равенство (после трёх умножений)
if(std::abs(d) < std::numeric_limits<float>::epsilon()*3)
{
x1re=x2re=-b/(2*a);
x1im=x2im=0;
}
else if(d<0)
{
x1re=x2re=-b/(2*a);
x1im=-(x2im=-std::sqrt(-d)/(2*a));
}
else
{
x1re=(-b-sqrt(d))/(2*a);
x2re=(-b+sqrt(d))/(2*a);
x1im=x2im=0;
}
}
// Шаблон для итерирования по количеству разрядов в дробной части
// значений с фиксированной точкой.
// В связи с невозможностью (запрещено стандартом) шаблонам функций иметь
// частичные специализации итерирации выполняются перегруженным operator()
// в шаблоне класса.
// В тесте не будем параметризировать тип параметров
template <__int64 a, __int64 b, __int64 c, int step>
class testOut
{
static const Qwur<__int64, a, b, c, step> qw;
public:
void operator()() const
{
std::cout<<tab<<qw.x1Re<<tab<<qw.x2Re<<tab<<qw.x1Im<<tab<<qw.x2Im<<std::endl;
testOut<a*10, b*10, c*10, step+1>()(); // Следующая точность
}
};
// Специализация теста для остановки рекурсии.
// Точность в семь знаков является по факту предельной для __int64.
// Однако и точность float неплоха для фиксированной точки в compile-time
template <__int64 a, __int64 b, __int64 c>
struct testOut<a, b, c, 7>
{
void operator()() const {}
};
int main()
{
float x1re, x1im, x2re, x2im;
testOut<7, 70, 175, 0>()();// Нулевой дискриминант
qwur (7, 70, 175, x1re, x1im, x2re, x2im);
std::cout<<tab<<x1re<<tab<<x2re<<tab<<x1im<<tab<<x2im<<std::endl<<std::endl;
testOut<7, 28,-315, 0>()();// Положительный дискриминант
qwur (7, 28, -315, x1re, x1im, x2re, x2im);
std::cout<<tab<<x1re<<tab<<x2re<<tab<<x1im<<tab<<x2im<<std::endl<<std::endl;
testOut<4, 16, 180, 0>()();// Отрицательный дискриминант
qwur (4, 16, 180, x1re, x1im, x2re, x2im);
std::cout<<tab<<x1re<<tab<<x2re<<tab<<x1im<<tab<<x2im<<std::endl<<std::endl;
testOut<1, 1, -1, 0>()(); // "Золотое сечение" (под корнем не полный квадрат)
qwur (1, 1, -1, x1re, x1im, x2re, x2im);
std::cout<<tab<<x1re<<tab<<x2re<<tab<<x1im<<tab<<x2im<<std::endl<<std::endl;
/*
К слову сказать, без реализации математического округления результат
последнего выводимого compile-time const отличался от run-time computable
*/
} |
Всего записей: 613 | Зарегистр. 08-08-2006 | Отправлено: 13:26 12-01-2007 | Исправлено: Qraizer, 13:49 12-01-2007 |
|
