Арифметика нового поколения: реализация числового формата Posit. Часть II

В первой части статьи мы разобрали что такое Posit, зачем он нужен и определились с его структурой. В этой части мы, как и обещали, расскажем о нашей реализации Posit на D.

Для реализации типа Posit!(NBITS, ES) нам потребуется некоторая предварительная подготовка: требуется разместить в файле модуля шаблон-генератор универсального свойства и описанные ранее битовые операции, которые будут использоваться далее постоянно. После этого, для того, чтобы иметь возможность строить Posit из double и приводимых к нему типов, а также для того, чтобы иметь полный доступ к структурным элементам Posit (знак, режим, экспонента и дробная часть) потребуется определить две вспомогательные структуры данных. Эти структуры данных будут использоваться только для внутреннего использования внутри библиотеки, реализующей тип Posit, и не предназначены для использования пользователями библиотеки.

DoubleComponents

Первая такая структура называется DoubleComponents и предназначена она для извлечения структурных элементов double, а именно, полей ответственных за знак (Sign), экспоненту (Exponent) и дробную часть (Fraction). Управление доступом к этим полям осуществляется с помощью автоматически сгенерированных сеттеров и геттеров, а самая важная часть — извлечение данных полей осуществляется в методе fromDouble, который принимает некое действительное число и с помощью ряда интересных битовых операций, выполняет выделение интересующих нас компонентов.

Структура данных DoubleComponents описана следующим образом:

class DoubleComponents
{
	mixin(addProperty!(long, "Value"));
	mixin(addProperty!(long, "Sign"));
	mixin(addProperty!(long, "Exponent"));
	mixin(addProperty!(long, "Fraction"));
		
	this()
	{
			
	}
		
	// extract data from double
	auto fromDouble(double number)
	{
		// bitwise conversion to ulong
		_value = *(cast(ulong*) &number);
		// sign extraction
		_sign = _value >> 63;
		// exponent extraction
		_exponent = ((_value & ((1UL << 63) - 1)) >> 52) - 1023;
		// fraction extraction
		_fraction = (1UL << 52) | (_value & ((1UL << 52) - 1));
	}
}

Выделение компонентов начинается с «грязного хака», основанного на прямой работе с памятью: число в формате double, так же как и Posit, представляет собой целочисленный тип-контейнер, в котором присутствуют определенные битовые поля. Структура этих полей описана в стандарте IEEE 754, с которым вы можете ознакомится по ссылке, мы же расскажем интересный момент, связанный с прямым доступом к числовому контейнеру.

Как говорилось ранее, double — это на самом деле, некое целое число, размерностью 64 бита, которое интерпретируется компьютером, как число с плавающей точкой. Но, ни одна операция приведения типа в D не поможет получить это целочисленное представление: и cast, и to (и даже специализированный union из стандартной библиотеки для дробных чисел) всего лишь обрежут дробное число до подходящего целого. И потому, мы решили пойти на хитрость: сначала мы получим адрес памяти по которому хранится искомый double, а затем с помощью приведения к указателю целочисленного типа, получим ссылку на якобы целое число. Дальнейшее просто: теперь достаточно разыменовать указатель типа ulong*, чтобы получить нужное нам целочисленное значение.

Следующим шагом выступает декодирование полей double, схему расположения которых вы можете увидеть на схеме:

При декодировании первым извлекается знак числа, для чего осуществляется доступ к биту под номером 63, схема интерпретации которого как знака точно такая же как в Posit. Это сделано тривиальной операцией сдвига.

После определения знака требуется вычленить саму экспоненту (или порядок числа), которая в double следует сразу за знаком и занимает 11 битов: именно для этого применяется целый комплекс битовых операций по удалению знака и итоговым сдвигам на 52 бита вправо. Помимо самих битовых операций, получение экспоненты привлекает в выражение вычитание числа 1023, поскольку согласно формату double, с его помощью описывается порядок и от этого числа идет его отсчет.

Извлечение дробной части (или мантиссы) осуществляется похожим образом, но занимает она ровно 52 бита, но простого сдвига на 52 позиции влево недостаточно, поскольку сама мантисса хранится в виде числа, деленного на 252 с прибавленной единицей. Поэтому требуется еще один такой же сдвиг, для декодирования величины мантиссы.

PositComponents

Следующая структура данных обеспечит библиотеку информацией о внутренней структуре Posit и называется она PositComponents. В нее с помощью нашего генератора свойств добавляются сеттеры и геттеры: для начального бита режима (RegimeSign), знака числа (Sign), самого режима (Regime), экспоненты (Exponent), дробной части (Fraction), а также для размера самих компонентов в битах (RegimeLength, ExponentLength и FractionLength соответственно). Также в PositComponents входит метод fromLong, который принимает целое число (тот самый контейнер под Posit) и интерпретируя его как битовое поле, разбирает на соответствующие компоненты. Каким образом это происходит, вы можете разобраться использовав описание самого формата Posit данного нами.

Структура данных PositComponents описана следующим кодом:

// representation of internal posit structure
class PositComponents
{
	// 0 or 1 in begin of regime
	mixin(addProperty!(long, "RegimeSign"));
	// sign bit
	mixin(addProperty!(long, "Sign"));
	// regime
	mixin(addProperty!(long, "Regime"));
	// exponent
	mixin(addProperty!(long, "Exponent"));
	// fraction
	mixin(addProperty!(long, "Fraction"));
	// regime length
	mixin(addProperty!(long, "RegimeLength"));
	// exponent length
	mixin(addProperty!(long, "ExponentLength"));
	// fraction length
	mixin(addProperty!(long, "FractionLength"));
		
	this(){}
		
	// extract component from long value (value is bit pattern)
	auto fromLong(long numberOfBits, long exponentSize)(long number)
	{
		long _number = number;
		_sign = getBit(_number, numberOfBits  - 1);
			
		if (_sign == 1)
		{
			_number = twoComplement(_number, numberOfBits);
		}
			
		_regimesign = getBit(_number, numberOfBits - 2);
						
		if (_regimesign == 0)
		{
			_regimelength = numberOfBits - lastSettedBit(_number) - 1;
		}
		else
		{
			_regimelength = numberOfBits - lastUnsettedBit(_number) - 1;
		}	
			
		_exponentlength = max(0, min(exponentSize, numberOfBits - 1 - _regimelength));
		_fractionlength = max(0, numberOfBits - 1 - _regimelength - _exponentlength);
			
		if (_regimesign == 0)
		{
			_regime = - _regimelength + 1;
		}  
		else
		{ 
			_regime = _regimelength - 2;
		}
			
		_exponent = extractBitField(_number, _exponentlength, _fractionlength) << (exponentSize - _exponentlength);
		_fraction = removeTrailingZeroes(
			setBit(
				extractBitField(_number, _fractionlength, 0), _fractionlength)
			);
	}
}

Эта структура имеет ключевое значение в библиотеке, так как активно используется для формирования Posit-контейнеров.

PositContainer

Сам контейнер для Posit является целочисленным типом, но с переменным размером, поскольку само определение формата говорит о том, что Posit может быть разного размера. Чтобы учесть это, потребуется небольшой шаблон, который в зависимости от заданного размера (в битах) сформирует нужный тип под контейнер и выглядит шаблон так:

// posit container
template PositContainer(uint size)
{
	static if (size == 8)
		alias PositContainer = ubyte;
	else static if (size == 16)
		alias PositContainer = ushort;
	else static if (size == 32)
		alias PositContainer = uint;
	else static if (size == 64)
		alias PositContainer = ulong;
	else
		static assert(0);
}

Здесь используется конструкция static if, которая позволяет на этапе компиляции проверить некоторые условия и в зависимости от них сформировать некоторый код. В нашем случае, мы создаем псевдоним PositContainer, который будет указывать на целочисленный тип, в котором будет хранится битовая строка Posit. Также, на этом этапе решено было избавиться от нестандартных размеров Posit и ограничиться только встречающимися в открытой документации на формат.

Posit!(NBITS,ES)

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

Определяем класс для Posit, в который добавляем несколько закрытых полей, которые (за исключением поля components типа PositComponents) будут закрыты для изменения и сформированы уже на этапе компиляции. Этими полями будут: Unum — псевдоним для типа контейнера (результат уже известен из шаблона), свойства get/setUnum для доступа к нему, NBITS — размер контейнера в битах, ES — максимальный размер экспоненты в битах, USEED — масштабирующий фактор useed, MINPOS — битовый паттерн минимально возможного положительного значения, MAXPOS — битовый паттерн максимально возможного положительного значения, NOT_A_REAL — битовый паттерн для обозначения бесконечности, ZERO — битовый паттерн для обозначения нуля, NPAT — общее количество всех битовых паттернов. Не все из этих полей на момент написания этой статьи реально использованы в библиотеке, но большая часть из них действительно применяется. Часть полей планировалось использовать для обработки некоторой арифметики, но к сожалению, описание ее реализации не было найдено, а неиспользуемые поля остались на будущее — для возможного расширения функционала.

Полный код класса Posit!(NBITS,ES) не публикуем в статье по причине его большого размера в строках кода, можете ознакомиться с ним в репозитории.

Мы не будем раскрывать подробности того, как все это работает, поскольку все это достаточно сложно и сопряжено с математикой. Также, описывать работу этого кода бессмысленно в силу того, что мы так и не нашли официальной документации по математическим операциям с Posit и в роли отправной точки использовалась библиотека PySigmoid. К этой библиотеке у нас возник ряд вопросов, которые, увы, остались без ответа, а потому наш код можно рассматривать как некий минимальный порт с Python, доведенный до ума в D. Несмотря на это, мы можем рассказать про присутствующие в классе методы и их функционал, в особенности тот, что пришлось сделать самим.

Ядром всего класса являются по сути несколько методов: construct — метод, формирующий данные в контенере из величины знака (sign), величины числа (scale) и дробной части (fraction); decode — метод, выделяющий компоненты Posit и сохраняющий их в структуру PositComponents; fromBits — метод, записывающий в контейнер некоторую строку бит (т.е. битовый паттерн всего Posit); fromDouble — метод, создающий число Posit из double.

Арифметические операции для типа реализованы были посредством стандартного способа перегрузки операций в D и рассчитаны на использование только с Posit-типом (размеры числа и максимальные размеры экспонент должны совпадать!), и были портированы из упоминавшейся библиотеки с некоторой доработкой: в исходной библиотеке не были учтены некоторые условия для умножения/деления, что могло бы привести к печальным результатам. При этом поддерживаются: унарный минус, сложение, вычитание, деление, нахождение обратной величины.

Операции сравнения для Posit реализованы через метод toUnsignedInteger и перегрузку opCmp: дело в том, что величины типа Posit легко и просто сравниваются, для этого всего лишь нужно сравнить между собой целочисленные контейнеры их хранящие, предварительно переведя контейнеры в формат целочисленного числа со знаком. Вот эту часть нам пришлось сделать самим, поскольку вариант из Python не будет работать в D: D — статически типизированный язык, и поэтому пришлось разработать шаблон SignedType, который по беззнаковому типу генерирует соответствующий ему тип со знаком. Данный шаблон используется в операции конверсии, без нее сравнение двух чисел Posit не приведет к успеху. Реализация перегрузки opCmp сделана нами крайне своеобразно: данная операция предполагает, что будут сравниваться между собой объекты (класс Object, корневой для всех объектов в D), однако, мы этого не планировали. Из-за этого, пришлось идти следующим путем: сравниваем типы исходного числа и поступившего с помощью typeid, вводим явное преобразование полученного объекта в Posit!(NBITS, ES), обе величины сравнения переводим в целое со знаком, а после этого всего сравниваем два числа со знаком и возвращаем результат в виде целого числа (1 — первое число больше; 0 — числа равны; -1 — первое число меньше). Таким образом, сравнение работает только между Posit, хотя при необходимости это можно изменить.

Для работы со сравнениями можно использовать не только реализованные в перегрузке <, >, =; но и ряд методов, которые мы реализовали самостоятельно: isMin — является ли минимальным положительным значением; isMax — является ли максимальным положительным значением; isNaR — является ли бесконечностью; isValid — является ли заданное значение корректным с точки зрения формата (особенно, что касается той битовой строки, из которой вы захотели создать Posit, в обход методов fromDouble/fromInteger); isZero — является ли нулем. Эти методы реализованы нами в полном объеме и напрямую использованы в реализации арифметики.

Вишенкой в нашем порте является то, чего не было реализовано в PySigmoid, но что определенно упростило отладку кода. Я говорю, про тот метод который помещен в секцию version(linux) — toBinaryFormatted. Этот метод сработает только на Linux системах, и вот почему. Дело в том, что ряд библиотек для работы с Posit, к примеру, SoftPosit предлагают очень интересную вещь — в них с помощью похожего по наименованию метода можно вывести битовую строку Posit с красивым форматированием, о котором я упоминал в первой части. Такое форматирование окрашивает цифры в определенные цвета: биты знака окрасятся красным, биты режима — желтым, биты экспоненты — синим, а биты дробной части — черным. А работает данное чудо на escape-последовательностях[здесь ссылка на статью по окраске терминала из блога] эмулятора терминала, что работает в Linux и некоторых других Unix-подобных системах, но не работает в Windows (но можно использовать эмуляторы линуксового терминала, например Git Bash).

Поскольку данная особенность доступна только на некоторых системах, а вывод Posit чисел в PySigmoid нас не устраивал (причина проста: не видно, какого типа контейнер, что в нем, и какова размерность; там Posit отображается так, как будто это float или double), не хочется оставаться без информации о Posit, то мы перегрузили метод toString так, что он стал выводить Posit в виде строки следующего формата:

Posit{размервбитах,максимальныйразмерэкспонентывбитах}(значение в контейнере в виде шестнадцатеричного числа)[эквивалентное значение в double]

В довершение всего вышесказанного, мы закрыли класс Posit!(NBITS,ES), сделав его приватным внутри модуля, и оставили псевдонимы для самых используемых типов: Posit8, Posit16, Posit32, Posit64 и SoftFloat (это экспериментальный тип Posit8 с двумя битами экспоненты), в итоге получается, что пользователю извне доступны только эти типы, общее надмножество остается недоступным.

Со всеми нюансами вы можете ознакомится в репозитории на GitHub.
Для демонстрации воспользуемся примером из этого репозитория:

import std.stdio;

import aftermath;

// SoftFloat is alias for  Posit!(8, 2): 8-bit Posit with 0 bits of exponent
auto factorial(SoftFloat p)
{
	SoftFloat n = new SoftFloat;
	n.fromInteger(1);
	
	if (p < n)
	{
		return n;
	}
	else
	{
		return p * factorial(p - n);
	}
}

// simple implementation of sinus function based on Taylor-McClaren series
auto sinus(SoftFloat p)
{
	SoftFloat r = new SoftFloat;
	
	for (int i = 0; i < 9; i++)
	{	
		SoftFloat o = new SoftFloat;
		o.fromInteger((-1) ^^ i);
		
		SoftFloat x = new SoftFloat;
		x.fromDouble(1.0);
		
		for (int j = 0; j < (2 * i + 1); j++)
		{
			x = x * p;
		}
		
		SoftFloat n = new SoftFloat;
		n.fromInteger(2 * i + 1); 
		
		//writeln("i ", i, " ", r);
		r = r + (o * (x / factorial(n)));
	}
	
	return r;
}

// fast sigmoid function
auto sigmoid(SoftFloat p)
{
	SoftFloat pp = new SoftFloat;
	
	auto n = p.toUnsignedInteger;
	n = n ^ (1 << (n.sizeof * 8 - 1));
	n >>= 2;
	
	pp.fromBits(n);
	return pp;
}




void main()
{
	// posit initialization
	SoftFloat p = new SoftFloat;
	// set posit to zero
	p.fromDouble(0);
	// writeln in string form
	p.writeln;
	// sigmoid result
        p.sigmoid.writeln;
	// set value of posit to 0.5
	p.fromDouble(0.5);
	// writeln new posit value
	p.writeln;
	// writeln in nice colored form (only on Linux)
	p.toBinaryFormatted.writeln;
	// get sine of 0.5 in posit
	p.sinus.writeln;
	// nice color form
	//p.sinus.toBinaryFormatted.writeln;
}

Что здесь происходит?

Мы используем SoftFloat, т.е. экспериментальный Posit тип для некоторых вычислений. Сначала мы определяем функцию факториала, которая описывается в простом рекурсивном стиле, и в этом нет ничего особенного: просто очередной новый числовой тип. Поскольку в библиотеке aftermath (так мы решили назвать наш порт) реализованы все базовые арифметические операции, то решено было попробовать нечто новое: реализовать на пробу какую-нибудь тригонометрическую функцию. В качестве такой функции был выбран синус, поскольку мы очень хорошо помним его разложение в ряд Макклорена. В разложении синуса в такой ряд используется факториал и возведение в степень, но последняя операция не была нами реализована. Поэтому есть дополнительный цикл в реализации рядов, который просто выполняет постоянное умножение, так как степень тут целочисленная. Поскольку, факториал уже реализован, то никаких остальных операций реализовать дополнительно не нужно: остается лишь последовательно вычислять члены ряда и суммировать их в некоторый аккумулятор.

Данный пример показывает, что все необходимое для вычислений у нас уже имеется, а все остальное реализуется из базовой арифметики. Тем не менее, в примере есть еще один момент: в оригинальной статье про сравнение Posit и float упоминалось, что при использовании формата float некоторые операции могут быть довольно накладными, а в качестве примера упоминалась сигмоидная функция, которую очень любят в нейронных сетях. Ради интереса мы решили создать такую функцию прямо по описанию: берем число Posit переключаем у него последний бит и сдвигаем на 2 вправо. Именно это и делается в реализации функции sigmoid.
Дальнейшее просто: мы создаем новое число и проделываем с ним те операции, которые смогли реализовать с одновременным выводом результатов на печать:

На этом мы заканчиваем эту непростую статью и просим тех читателей, кто дошел до этого момента и кто имеет желание/понимание, обратить внимание на этот сложный проект и принять горячее участие в разработке мощной библиотеки для Posit в D. Также, хотелось бы иметь Posit в стандартной библиотеке языка, чтобы быть на стороне прогресса…

Принять прямое участие в обсуждении наших проектов, а также ваших идей можно а нашем Discord: https://discord.gg/Xxyw4xW

Добавить комментарий