Введение

Библиотека AMP (ALGLIB Muliple Precision) - это объектно-ориентированный интерфейс к библиотеке арифметики высокой точности MPFR, созданный в рамках проекта ALGLIB (см. www.alglib.net и alglib.sources.ru). Для объяснения причин, по которым потребовалось создавать ещё один интерфейс к MPFR, следует немного рассказать о проекте ALGLIB.

Проект ALGLIB заключается в создании библиотеки численного анализа, в которой каждый алгоритм доступен в виде программ на нескольких языках программирования. В проекте ALGLIB эта задача решена за счет использования автоматического перевода. Алгоритмы хранятся в виде псевдокода, который переводится при помощи автоматического транслятора на требуемый язык. Синтаксис псевдокода позволяет осуществлять автоматический перевод на наиболее популярные языки программирования (C++, C#, Pascal и другие). Особый интерес представляет то, что в ряде случаев программа на псевдокоде может быть автоматически модифицирована таким образом, что она сможет использовать в своей работе арифметику высокой точности. Библиотека AMP при этом используется в качестве интерфейса к MPFR, заменяющего собой стандартный тип данных double.

То, что программы, использующие AMP, изначально были написаны для работы с вещественной арифметикой 64-битной точности, повлияло на дизайн библиотеки. Этим были обусловлены такие особенности, как определение точности на этапе компиляции (в AMP точность задается при помощи параметризации класса-шаблона) и запрет на операции со смешанной точностью (например, сложение 256-битного и 512-битного чисел). Эти ограничения функциональности необходимы для того, чтобы библиотеку было можно использовать в проекте ALGLIB.

Следует отметить, что хотя библиотека AMP в первую очередь предназначена для проекта ALGLIB, она может представлять интерес и в качестве "просто интерфейса к MPFR". В частности, автоматизированная сборка мусора может оказаться полезной в ряде приложений. В будущем планируется развивать библиотеку именно в этом направлении, т.к. относящаяся непосредственно к ALGLIB часть библиотеки уже достаточно проработана и не нуждается в дальнейшем расширении.

Лицензия

Библиотека AMP доступна на условиях тройной лицензии MPL 1.1 / GPL 2.0 / LGPL 2.1

Разъяснение. Предоставление библиотеки на условиях тройной лицензии обозначает, что пользователь может использовать AMP на условиях любой из трех лицензий по своему выбору. При распространении AMP и/или её модификаций пользователю доступны два варианта поведения. Во-первых, он может распространять её на тех же условиях тройной лицензии, оставляя за получателями библиотеки право использовать любую из трех разрешенных лицензий по своему усмотрению. Во-вторых, он имеет право распространять библиотеку на условиях лицензии GPL 2.0 (или, по своему выбору, LGPL 2.1), запретив распространение на условиях других лицензий.

Эти условия распространения позволяют использовать библиотеку AMP как в проприетарных или просто несвободных программах (что допускается условиями лицензии MPL), так и в ПО, распространяющемся по лицензии GPL/LGPL. По мнению автора, это способствует максимально широкому распространению и использованию библиотеки.

Подключение и использование

Для использования библиотеки AMP необходимо:

Основные концепции

Вещественные числа. MPFR (и, следовательно, AMP) практически полностью соответствует стандарту IEE-754, содержащему требования к реализации вещественной арифметики. Пользователю доступны следующие типы чисел: положительный и отрицательный нули (signed zeros), конечные числа, положительная и отрицательная бесконечности, NaN. Денормализованные числа не поддерживаются.

Операции с числами и округление. Библиотека MPFR поддерживает четыре типа округления: округление вверх, вниз, к нолю и к ближайшему представимому значению. Библиотека AMP всегда использует режим округления к ближайшему представимому значению. Таким образом, все операции с вещественными числами проводятся с использованием точного округления (exact rounding). Появление в ходе вычислений специальных значений (бесконечностей и NaN) не приводит к генерации исключения - вычисления продолжаются в нормальном режиме. Операции с бесконечностями и нулями со знаком осуществляются по правилам арифметики пределов.

Точность. Точность вещественного числа - это число бит в мантиссе, включая старший (ненулевой бит). Например, если установить точность вещественного числа равной 53, то получится тип данных, имитирующий стандартный тип double. Библиотека AMP позволяет устанавливать в качестве точности любое число, не меньшее 32 (размер типа данных signed long), включая и числа, не являющиеся степенями двойки. Это ограничение введено из соображений совместимости - многие математические программы полагаются на то, что целое значение может быть присвоено вещественному без потери точности.

Определение точности на этапе компиляции. Особенностью библиотеки AMP является то, что класс, позволяющий осуществлять операции с вещественными числами высокой точности, является классом-шаблоном, параметром которого является целое число, задающее точность хранящегося значения. Таким образом, точность, с которой осуществляются вычисления, становится известна ещё на этапе компиляции - при параметризации шаблона, и не может быть изменена в ходе выполнения программы. Это существенное ограничение функциональности библиотеки было введено сознательно. Причиной является то, что основное назначение библиотеки AMP - применение в рамках проекта ALGLIB. В проекте ALGLIB библиотека AMP используется программами, являющимися результатом работы автоматического транслятора. Генерируемый транслятором код не способен работать с вещественными числами, чья точность меняется динамически, поскольку при этом возникает ряд проблем (например, какую точность следует использовать для промежуточных результатов, если одни элементы матрицы заданы с точностью 128 бит, а другие - с точностью 256 бит). По этой причине было решено ввести вышеуказанное ограничение.

Операнды смешанной точности. С++ предоставляет несколько вещественных типов: float, double, long double. При этом, если аргументами какой-либо операции являются два числа различного типа, то обычно осуществляется неявное преобразование от меньшего типа к большему (так, сумма float и double имеет тип double). Библиотека AMP не позволяет осуществлять неявное приведение типов при операциях с операндами смешанной точности. Например, сложение amp::ampf<128> и amp::ampf<144> вызовет ошибку компиляции - потребуется в явной форме привести один из типов к типу другого операнда. Отчасти это продиктовано ограничениями С++ (которые в принципе можно обойти), отчасти - убеждением разработчика в том, что в практике численного анализа потребность в подобной функциональности возникает очень редко.

Сборка мусора. Поскольку в MPFR память под вещественные числа выделяется динамически, инициализация нового вещественного числа занимает сравнительно много времени. Для ускорения работы библиотека AMP не освобождает использованные переменные, а помещает их в список, из которого осуществляется выборка, когда требуется создать новое вещественное число.

Отсутствие многопоточности. В настоящее время библиотека AMP не поддерживает многопоточность. Библиотека может работать в многопоточной среде, но не может использоваться двумя и более потоками одновременно.

Классы исключений

Библиотека AMP включает в себя следующие классы исключений:

  1. amp::incorrectPrecision - генерируется при попытке инициализировать число с указанием недопустимой точности (менее 32 или больше ограничений, устанавливаемых библиотекой MPFR)
  2. amp::invalidConversion - генерируется при попытке привести к целому виду при помощи функций trunc/round/ceil/floor число, не помещающееся в разрядную сетку типа данных signed long
  3. amp::internalError - генерируется при внутреннем сбое в библиотеке AMP (в нормальных условиях это не должно происходить)

Классы amp::mpfr_record и amp::mpfr_storage

Эти классы - служебные классы, используемые библиотекой для хранения указателя на вещественное число и как интерфейс менеджера памяти. Хотя эти классы доступны программисту, он не должен использовать их, поскольку они не предназначены для использования за пределами библиотеки AMP. Правила работы с этими классами не документированы, а несоблюдение этих правил способно разрушить систему выделения памяти или нарушить нормальный процесс вычислений.

Класс amp::ampf и операции с ним

Общие принципы

Класс amp::ampf является классом-шаблоном, инкапсулирующим вещественное число в формате библиотеки MPFR. Класс полностью берет на себя такие функции, как выделение памяти для вещественного числа, освобождение памяти и уборка мусора. Для объектов этого типа определены арифметические операции и ряд функций, копирующих функциональность стандартной библиотеки языка С. В большинстве случаев, код, использующий класс amp::ampf, почти не отличается от кода, работающего с вещественными числами.

Члены класса amp::ampf


    template<unsigned int Precision>
    class ampf
    {
    public:
        ~ampf();

        ampf ();        
        ampf (long double v);
        ampf (double v);
        ampf (float v);
        ampf (signed long v);
        ampf (unsigned long v);
        ampf (signed int v);
        ampf (unsigned int v);
        ampf (signed short v);
        ampf (unsigned short v);
        ampf (signed char v);
        ampf (unsigned char v);
        ampf (const std::string &s);
        ampf (const char *s);
        
        ampf(const ampf& r);
        template<unsigned int Precision2>
        ampf(const ampf<Precision2>& r);

        ampf& operator= (long double v);
        ampf& operator= (double v);
        ampf& operator= (float v);
        ampf& operator= (signed long v);
        ampf& operator= (unsigned long v);
        ampf& operator= (signed int v);
        ampf& operator= (unsigned int v);
        ampf& operator= (signed short v);
        ampf& operator= (unsigned short v);
        ampf& operator= (signed char v);
        ampf& operator= (unsigned char v);
        ampf& operator= (const char *s);
        ampf& operator= (const std::string &s);
        ampf& operator= (const ampf& r);
        template<unsigned int Precision2>
        ampf& operator= (const ampf<Precision2>& r);
        
        template<class T>
        ampf& operator+=(const T& v);
        template<class T>
        ampf& operator-=(const T& v);
        template<class T>
        ampf& operator*=(const T& v);
        template<class T>
        ampf& operator/=(const T& v);
        
        mpfr_srcptr getReadPtr() const;
        mpfr_ptr getWritePtr();
        
        bool isFiniteNumber() const;
        bool isPositiveNumber() const;
        bool isZero() const;
        bool isNegativeNumber() const;
        const ampf getUlpOf();

        double toDouble() const;
        std::string toHex() const;
        std::string toDec() const;
        
        static const ampf getUlpOf(const ampf &x);
        static const ampf getUlp();
        static const ampf getUlp256();
        static const ampf getUlp512();
        static const ampf getMaxNumber();
        static const ampf getMinNumber();
        static const ampf getAlgoPascalEpsilon();
        static const ampf getAlgoPascalMaxNumber();
        static const ampf getAlgoPascalMinNumber();
        static const ampf getRandom();
    };

Объявление переменных и инициализация

Класс amp::ampf является классом-шаблоном, параметризируемым по точности хранящегося в нем вещественного значения. Таким образом, перед использованием этого класса необходимо определить, какая именно точность требуется. Для инициализации класса определен ряд конструкторов, принимающих как объекты того же класса, так и целые числа, вещественные числа стандартных типов (float, double, long double) и строки, в которых число записано в дробной (например, "12.34") или экспоненциальной форме, при этом мантисса может быть как в десятичной системе исчисления, так и в шестнадцатеричной системе (см. описание библиотеки MPFR).

Примеры:


    int i = 3;
    double d = 2.5;
    amp::ampf<128> a(i), b(d), c("2.33");
    a = "1.56E-3";
    b = a;

Преобразование в текстовую форму

Для преобразования числа в текстовую форму служат функции-члены toDec() и toHex(). Они возвращают экспоненциальное представление числа в десятичной и шестнадцатеричной системе исчисления (при этом шестнадцатеричное преобразование осуществляется быстрее, чем десятичное, что делает его предпочтительным в тех случаях, когда требуется повышенное быстродействие). Какие-либо возможности по форматированию результата отсутствуют: число выдается с максимальным количеством значащих цифр.

Пример:


    amp::ampf<128> a;
    a = "11";
    printf("%s\n", amp::sqrt(a).toDec().c_str());

Проверка свойств

Функции-члены isFiniteNumber(), isPositiveNumber(), isZero(), isNegativeNumber() позволяют проверить, является ли число конечным значением (в противоположность одниму из специальных значений), положительным числом, нулем или отрицательным числом.

Информационные функции

Ряд функций предназначен для получения информации о числовых свойствах того или иного типа данных:

Доступ к указателю mpfr_t

Доступ к указателю mpfr_t предоставляют две функции: getReadPtr() и getWritePtr(). Функция getReadPtr() позволяет получить указатель только для чтения, который можно передавать в функции библиотеки MPFR в качестве аргумента. Функция getWritePtr() позволяет получить указатель для записи, который можно передавать в функции библиотеки MPFR в качестве места для сохранения результата вычислений.

Следует отметить, что библиотека AMP использует технику copy-on-write, то есть указатель, полученный методом getReadPtr() может быть общим для нескольких объектов. Таким образом, указатель, полученный методом getReadPtr() ни при каких условиях не должен использоваться в операциях, которые могут модифицировать его содержимое.

Указатель, полученный методом getWritePtr(), может быть использован в операциях, меняющих его содержимое, но эти операции не должны менять его точность. Следует учитывать, что из-за использования техники copy-on-write, указатель, полученный при вызове getWritePtr(), может отличаться от указателя, полученного ранее при вызове getReadPtr().

Генерация случайных чисел

Для генерации случайных чисел служит статическая функция getRandom(). При первом вызове эта функция инициализирует генератор случайных чисел используя системное время.

Арифметические операции

Класс amp::ampf переопределяет операции сравнения, сложения, вычитания, умножения и деления. Результат округляется по правилам точного округления (exact rounding). Операндами могут быть два объекта amp::ampf одной и той же точности. В выражении допускается смешивать объекты типа amp::ampf и стандартные типы данных. Смешивать объекты типа amp::ampf с разной точностью можно, только если один из операндов явным образом приводится к типу другого.

Пример:


    amp::ampf<128> a(3), b(4), c(99);
    amp::ampf<256> d(1);

    a = b+c;          //
    b = c+2.3;        // допустимые выражения
    c = (2*a-5)/b;    //

    a = b+d;          // вызовет ошибку при компиляции -
                      // аргументы имеют разный тип

    a = b + amp::ampf<128>(d); // так правильно

    a = d;            // это тоже допустимо

Математические функции

Все приведенные в этом разделе функции являются членами пространства имен amp. Таким образом, при вызове, например, функции sqr, следует писать amp::sqr. Если явно не указано иное, то результат округляется по правилам точного округления.

template<unsigned int Precision>
const ampf<Precision> sqr(const ampf<Precision> &x)
Эта функция служит для возведения числа в квадрат.

template<unsigned int Precision>
const int sign(const ampf<Precision> &x)
Эта функция возвращает +1, если аргумент положителен, -1, если он отрицателен, 0, если он равен нолю.

template<unsigned int Precision>
const ampf<Precision> abs(const ampf<Precision> &x)
Эта функция возвращает модуль числа.

template<unsigned int Precision>
const ampf<Precision> maximum(const ampf<Precision> &x, const ampf<Precision> &y)
Эта функция возвращает большее из двух чисел.

template<unsigned int Precision>
const ampf<Precision> minimum(const ampf<Precision> &x, const ampf<Precision> &y)
Эта функция возвращает меньшее из двух чисел.

template<unsigned int Precision>
const ampf<Precision> sqrt(const ampf<Precision> &x)
Эта функция возвращает квадратный корень.

template<unsigned int Precision>
const signed long trunc(const ampf<Precision> &x)
Эта функция округляет число в сторону ноля и возвращает значение типа signed long. Если результат округления не помещается в разрядную сетку, то генерируется исключение amp::invalidConversion.

template<unsigned int Precision>
const ampf<Precision> frac(const ampf<Precision> &x)
Эта функция возвращает дробную часть числа. Дробная часть определяется, как x-trunc(x).

template<unsigned int Precision>
const signed long floor(const ampf<Precision> &x)
Эта функция округляет число вниз (в сторону минус бесконечности) и возвращает значение типа signed long. Если результат округления не помещается в разрядную сетку, то генерируется исключение amp::invalidConversion.

template<unsigned int Precision>
const signed long ceil(const ampf<Precision> &x)
Эта функция округляет число вверх (в сторону плюс бесконечности) и возвращает значение типа signed long. Если результат округления не помещается в разрядную сетку, то генерируется исключение amp::invalidConversion.

template<unsigned int Precision>
const signed long round(const ampf<Precision> &x)
Эта функция округляет число к ближайшему целому и возвращает значение типа signed long. Аргумент, располагающийся точно на границе, округляется в сторону бесконечности. Если результат округления не помещается в разрядную сетку, то генерируется исключение amp::invalidConversion.

template<unsigned int Precision>
const ampf<Precision> frexp2(const ampf<Precision> &x, mp_exp_t *exponent)
Эта функция разделяет число на мантиссу и экспоненту по основанию два. Если аргумент равен нолю, то мантисса и экспонента считаются равными нолю. Если аргумент не является конечным числом, то генерируется исключение amp::invalidConversion.

template<unsigned int Precision>
const ampf<Precision> ldexp2(const ampf<Precision> &x, mp_exp_t exponent)
Эта функция умножает переданный аргумент на два в степени exponent.

template<unsigned int Precision>
const ampf<Precision> pi()
Эта функция возвращает константу Pi.

template<unsigned int Precision>
const ampf<Precision> halfpi()
Эта функция возвращает половину константы Pi.

template<unsigned int Precision>
const ampf<Precision> twopi()
Эта функция возвращает удвоенную константу Pi.

template<unsigned int Precision>
const ampf<Precision> sin(const ampf<Precision> &x)
Эта функция возвращает синус аргумента.

template<unsigned int Precision>
const ampf<Precision> cos(const ampf<Precision> &x)
Эта функция возвращает косинус аргумента.

template<unsigned int Precision>
const ampf<Precision> tan(const ampf<Precision> &x)
Эта функция возвращает тангенс аргумента.

template<unsigned int Precision>
const ampf<Precision> asin(const ampf<Precision> &x)
Эта функция возвращает арксинус аргумента.

template<unsigned int Precision>
const ampf<Precision> acos(const ampf<Precision> &x)
Эта функция возвращает арккосинус аргумента.

template<unsigned int Precision>
const ampf<Precision> atan(const ampf<Precision> &x)
Эта функция возвращает арктангенс аргумента.

template<unsigned int Precision>
const ampf<Precision> atan2(const ampf<Precision> &y, const ampf<Precision> &x)
Эта функция возвращает арктангенс числа, равного отношению аргументов y/x. Функция возвращает корректный результат, даже если x равно 0.

template<unsigned int Precision>
const ampf<Precision> log(const ampf<Precision> &x)
Эта функция возвращает натуральный логарифм x.

template<unsigned int Precision>
const ampf<Precision> log2(const ampf<Precision> &x)
Эта функция возвращает двоичный логарифм x.

template<unsigned int Precision>
const ampf<Precision> log10(const ampf<Precision> &x)
Эта функция возвращает десятичный логарифм x.

template<unsigned int Precision>
const ampf<Precision> exp(const ampf<Precision> &_x)
Эта функция возвращает е в степени x.

template<unsigned int Precision>
const ampf<Precision> sinh(const ampf<Precision> &x)
Эта функция возвращает гиперболический синус x.

template<unsigned int Precision>
const ampf<Precision> cosh(const ampf<Precision> &_x)
Эта функция возвращает гиперболический косинус x.

template<unsigned int Precision>
const ampf<Precision> tanh(const ampf<Precision> &_x)
Эта функция возвращает гиперболический тангенс x.

template<unsigned int Precision>
const ampf<Precision> pow(const ampf<Precision> &x, const ampf<Precision> &y)
Эта функция возвращает x в степени y.

Векторы, матрицы и базовые операции линейной алгебры

Введение

Библиотека AP, входящая в состав проекта ALGLIB, включает в себя классы-шаблоны векторов и матриц, а также подпрограммы для осуществления над ними базовых операций линейной алгебры (аналог Level 1 BLAS). Поскольку библиотека AP реализует классы матриц, используя шаблоны, эти же классы можно использовать для работы с матрицами, состоящими из чисел высокой точности. По этой причине библиотека AMP не содержит классов матриц и векторов. Базовые операции линейной алгебры в библиотеке AP также реализованы с использованием шаблонов, но в этом случае библиотека AMP предлагает свой аналог, написанный специально для операций с числами высокой точности. Подпрограммы библиотеки AMP специально оптимизированы для снижения нагрузки на менеджер памяти, поэтому при осуществлении векторных операций над числами высокой точности следует использовать именно их.

Описание

template<unsigned int Precision>
ampf<Precision> vDotProduct(ap::const_raw_vector< ampf<Precision> > v1, ap::const_raw_vector< ampf<Precision> > v2)
Эта функция возвращает скалярное произведение двух векторов. Векторы задаются при помощи структур ap::const_raw_vector (эта структура определена в файле ap.h, более подробно см. описание библиотеки AP).

template<unsigned int Precision>
void vMove(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
Эта функция копирует вектор vSrc в вектор vDst.

template<unsigned int Precision>
void vMoveNeg(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
Эта функция копирует вектор vSrc в вектор vDst, умножая его на -1.

template<unsigned int Precision, class T2>
void vMove(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
Эта функция копирует вектор vSrc в вектор vDst, умножая его на константу alpha.

template<unsigned int Precision>
void vAdd(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
Эта функция добавляет к вектору vDst вектор vScr.

template<unsigned int Precision, class T2>
void vAdd(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
Эта функция добавляет к вектору vDst вектор vScr, умноженный на константу alpha.

template<unsigned int Precision>
void vSub(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
Эта функция вычитает из вектора vDst вектор vScr.

template<unsigned int Precision, class T2>
void vSub(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
Эта функция вычитает из вектора vDst вектор vScr, умноженный на константу alpha.

template<unsigned int Precision, class T2>
void vMul(ap::raw_vector< ampf<Precision> > vDst, T2 alpha)
Эта функция умножает вектор vDst на константу alpha.

Класс amp::campf и операции с ним

Общие принципы

Класс amp::campf реализует операции с комплексными числами высокой точности. Действительная и мнимая части хранятся в полях x и y, имеющих тип amp::ampf. С объектами этого типа можно осуществлять арифметические операции, операции сравнения, также для этого типа данных определены некоторые из стандартных математических функций.

Члены класса amp::ampf


    template<unsigned int Precision>
    class campf
    {
    public:
        campf():x(0),y(0);
        campf(long double v);
        campf(double v);
        campf(float v);
        campf(signed long v);
        campf(unsigned long v);
        campf(signed int v);
        campf(unsigned int v);
        campf(signed short v);
        campf(unsigned short v);
        campf(signed char v);
        campf(unsigned char v);
        campf(const ampf<Precision> &_x);
        campf(const ampf<Precision> &_x, const ampf<Precision> &_y);
        campf(const campf &z):x(z.x),y(z.y);
        template<unsigned int Prec2>
        campf(const campf<Prec2> &z);

        campf& operator= (long double v);
        campf& operator= (double v);
        campf& operator= (float v);
        campf& operator= (signed long v);
        campf& operator= (unsigned long v);
        campf& operator= (signed int v);
        campf& operator= (unsigned int v);
        campf& operator= (signed short v);
        campf& operator= (unsigned short v);
        campf& operator= (signed char v);
        campf& operator= (unsigned char v);
        campf& operator= (const char *s);
        campf& operator= (const std::string &s);
        campf& operator= (const campf& r);
        template<unsigned int Precision2>
        campf& operator= (const campf<Precision2>& r);

        ampf<Precision> x, y;
    };

Арифметические операции

Класс amp::campf переопределяет операции сравнения, сложения, вычитания, умножения и деления. Операндами могут быть два объекта amp::campf одной и той же точности. В выражении допускается смешивать объекты типа amp::campf, объекты типа amp::ampf и стандартные типы данных. Смешивать объекты типа amp::campf с разной точностью можно, только если один из операндов явным образом приводится к типу другого.

Стандартные функции

template<unsigned int Precision>
const ampf<Precision> abscomplex(const campf<Precision> &z)
Эта функция возвращает модуль комплексного числа.

template<unsigned int Precision>
const campf<Precision> conj(const campf<Precision> &z)
Эта функция возвращает число, сопряженное данному

template<unsigned int Precision>
const campf<Precision> csqr(const campf<Precision> &z)
Эта функция возвращает квадрат числа.