Числа с плавающей точкой (одинарной или двойной точности) печально известны за нелепые вычислительные ошибки. Например:

Вы можете увидеть больше примеров и даже оценок того, как часто возникают такие ошибки вычислений, здесь. В этой статье используется JavaScript, но написанное актуально для любого языка программирования, который реализует вещественные числа в формате с плавающей точкой, используя стандарт IEEE-754 (т.е. для любого современного языка программирования).

Почему это происходит?

Ресурсы компьютера ограничены. Например, нет способа сохранить число больше 255 или меньше 0 в одном байте. Если мы хотим сохранить знаковое значение, диапазон становится [-128; 127]. Вещественные числа в памяти компьютера также ограничены по точности. Целые числа вызывают переполнение, когда они выходят за пределы, но значения переменных типа с плавающей точкой вырождаются по мере увеличения их значений. Давайте рассмотрим пример:

import (
	"fmt"
)

func main() {
	fmt.Println(float32(16777215.0) == float32(16777216.0)) // false
	fmt.Println(float32(16777216.0) == float32(16777217.0)) // true
	fmt.Println(float32(16777216.0)) // 1.6777216e+07
	fmt.Println(float32(16777217.0)) // 1.6777216e+07
}

Как можно заметить, даже некоторые целые числа не могут быть адекватно представлены во float32, начиная с некоторой отметки. Тип float64 имеет более широкий диапазон целых чисел, вычисляемых без ошибок. Но что это за отметка? Ну, 16 777 216 = 2²³ × 2. А мантисса числа float32 имеет длину 23 бита. Действительно ли это та отметка, после которой значения целых чисел с плавающей точкой начинают вырождаться? Легко выяснить с помощью следующего кода:

for i := float32(0.0); i <= float32(16777216.0); i += float32(1.0) {
	if i == 16777215.0 {
		fmt.Printin("halt") // breakpoint here
	}

	if i == i+float32(1.0) {
		fmt.Println(i) // breakpoint here
	}
}

Этот небольшой блок инструкций ничего не печатает, пока значение переменной не достигнет 16 777 215. Когда она достигает 16 777 216, то сравнивается с 16 777 216 + 1, и это число не может быть представлено во float32. Поэтому код печатает значение счетчика. По той же причине увеличить его нельзя, и цикл будет продолжать печатать тот же вывод бесконечно.

То же самое можно воспроизвести и в отрицательном направлении. Вот как это выглядит:

for i := float32(0.0); i >= -float32(16777216.0); i -= float32(1.0) {
	if i == -16777215.0 {
		fmt.Printin("halt") // breakpoint here
	}
	
	if i == i-float32(1.0) {
		fmt.Println(i) // breakpoint here
	}
}

Мантисса числа float64 больше (52 бита в длину). Это значит, что вы можете свободно использовать целые числа с плавающей точкой, расположенные в диапазоне [-2⁵³; 2⁵³]. На самом деле, это довольно хороший диапазон, но да, то же самое правило применимо и здесь. Вот код:

package main

import (
	"fmt"
	"math"
)

func main() {
	fmt.Println(math.Pow(2, 53) == math.Pow(2, 53)-1) // false
	fmt.Println(math.Pow(2, 53) == math.Pow(2, 53)+1) // true
}

Почему это приемлемо?

- Добрый день. Я слышал, вы умеете быстро считать в уме. Давайте проверим. Итак, сколько будет 14 * 17?
- 230
- Но это неверно!
- Зато быстро.

Как с этим жить

Вы можете использовать числа с фиксированной точкой, тип Decimal или простое округление (если вы точно знаете, сколько цифр после точки должно быть именно в вашем случае). Если в вашем случае вам приходится иметь дело с денежными единицами, то имеет смысл применить шаблон Мартина Фаулера для валютных величин.

Если вам железно надо оставить float и вы не знаете, какова должна быть точность для округления, у меня есть для вас хитроумное решение.

Безопасные операции

Вычислительных ошибок не возникает при следующих операциях над двумя числами с плавающей точкой:

• сложение двух целых чисел в формате с плавающей точкой, значения которых ограничены вышеупомянутыми диапазонами ([-2²⁴, 2²⁴] для float32, [-2⁵³; 2⁵³] для float64); результат также должен быть в указанных пределах;

• вычитание двух целых чисел (те же правила);

• умножение двух целых чисел (те же правила);

• деление двух целых чисел (те же правила);

• умножение любого числа с плавающей точкой (целого или дробного) на 10 в степени n, где n ≥ 0, n — целое число;

• деление любого числа с плавающей точкой (целого или дробного) на 10 в степени n, где n ≥ 0, n — целое число;

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

Теория

Предположим, мы хотим прибавить 3,8 к 2,56. Если бы мы сделали это с помощью языка программирования, оператор вернул бы 6,359999999999999. Теперь попробуем сделать это, используя только безопасные операции. Для этого нам придется вспомнить немного базовой школьной математики:

Здесь мы просто складываем два целых числа и делим их на 10². То же самое применимо и к вычитанию. Однако умножение отличается, как можно видеть ниже:

Именно, все дело в нахождении длины дробной части числа и умножении ее на 10 в степени этой длины. Это позволяет работать с целыми числами вместо вещественных. В конце вычисления вам придется сделать результат дробным, используя деление на 10 в степени, которая зависит от оператора.

В случае сложения и вычитания степень — это максимальная длина дробной части среди обоих операндов. В умножении степень — это сумма длин дробных частей обоих операндов. Но с делением всё не так просто. Есть три случая:

• числитель имеет большую длину дробной части, чем знаменатель:
  

  

• знаменатель имеет большую длину дробной части, чем числитель:
  

  

• дробные части обоих чисел одинаковы по длине:
  

  

В первом случае мы делим результат. Во втором случае мы умножаем результат. В третьем случае мы оставляем его как есть.

Нахождение длины дробной части

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

import (
	"strconv"
	"strings"
)

func GetPartsLength(number float64) (int, int) {
	numStr := strconv.FormatFloat(number, 'f', -1, 64)
	numStr = strings.TrimPrefix(numStr, "-")
	parts := strings.Split(numStr, ".")

	if len(parts) > 1 {
		return len(parts[0]), len(parts[1])
	}

	if len(parts) == 1 {
		return len(parts[0]), 0
	}

	return -1, -1
}

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

func BenchmarkGetPartsLength(b *testing.B) {
	for i := 0; i < b.N; i++ {
		safecomp.GetPartsLength(-19367.5889)
	}
}

Результат не очень воодушевляющий:

Running tool: /usr/bin/go test -benchmem -run=A$ \
-bench ABenchmarkGetPartsLength$ ieee-754/safecomp

goos: linux
goarch: amd64
pkg: ieee-754/safecomp
cpu: AMD Ryzen 5 5600H with Radeon Graphics

BenchmarkGetPartsLength-12  7736520   263.0 ns/op  72 B/op   3 allocs/op

PASS ok ieee-754/safecomp 2.196s

Строковые операции, как правило, весьма затратны в плане использования ЦП и ОЗУ. Кроме того, они часто вызывают выделение динамической памяти. Должен быть другой способ, потому что strconv.FormatFloat должен каким-то образом знать, где разместить точку в строковом представлении числа типа float. Поэтому я решил посмотреть, что происходит внутри этой функции. Стэк вызовов выглядит следующим образом:

Функция ryuFtoaShortest реализует алгоритм Ryū. Она заполняет структуру типа decimalSlice:

type decimalslice struct {
	d      []byte
	nd, dp int
	neg    bool
}

nd означает «количество цифр». Оно хранит общую длину числа, исключая знак и точку. dp означает «десятичная часть». Оно хранит длину десятичной части, поэтому длина дробной части равна nd - dp.

Теперь пришло время заставить функции служить нашей цели. Они находятся в /usr/lib/go/src/strconv. Конечно, расположение зависит от пути установки, версии языка, ОС и т. д. Необходимые функции используют другие процедуры из пакета, поэтому было бы слишком утомительно копировать-вставлять только то, что нам нужно. Вот почему я решил перенести весь пакет в свой проект. И вот к чему это привело:

[zergon321@astroemeria ieee-754]$ go run main.go
package command-line-arguments
  imports ieee-754/alt
  alt/bytealg.go:10:8: use of internal package internal/bytealg not allowed

Упс. Что ж, я просто удалил файл bytealg.go и всё, что от него зависело. К счастью, функция genericFtoa к таковому не относилась. Все мои действия привели к написанию GetPartsLength, что показано ниже:

package alt

import "math"

func GetPartsLength(val float64, fmt byte, prec, bitSize int) (int, int) {
	var bits uint64
	var flt *floatInfo
	switch bitSize {
	case 32:
		bits = uint64(math.Float32bits(float32(val)))
		flt = &float32info
	case 64:
		bits = math.Float64bits(val)
		flt = &float64info
	default:
		panic("strconv: illegal AppendFloat/FormatFloat bitSize")
	}

	exp := int(bits>>flt.mantbits) & (1<<flt.expbits - 1)
	mant := bits & (uint64(1)<<flt.mantbits - 1)

	switch exp {
	case 1<<flt.expbits - 1:
		// Inf, NaN
		return -1, -1

	case 0:
		// denormalized
		exp++

	default:
		// add implicit top bit
		mant |= uint64(1) << flt.mantbits
	}
	exp += flt.bias

	var digs decimalSlice
	// Negative precision means "only as much as needed to be exact."
	shortest := prec < 0
	if shortest {
		// Use Ryu algorithm.
		var buf [32]byte
		digs.d = buf[:]
		ryuFtoaShortest(&digs, mant, exp-int(flt.mantbits), flt)
	} else if fmt != 'f' {
		// Fixed number of digits.
		digits := prec
		switch fmt {
		case 'e', 'E':
			digits++
		case 'g', 'G':
			if prec == 0 {
				prec = 1
			}
			digits = prec
		default:
			// Invalid mode.
			digits = 1
		}
		var buf [24]byte
		if bitSize == 32 && digits <= 9 {
			digs.d = buf[:]
			ryuFtoaFixed32(&digs, uint32(mant), exp-int(flt.mantbits), digits)
		} else if digits <= 18 {
			digs.d = buf[:]
			ryuFtoaFixed64(&digs, mant, exp-int(flt.mantbits), digits)
		}
	}

	return digs.dp, digs.nd - digs.dp
}

Это не идеально, но работает. Результат анализа производительности теперь намного лучше:

Running tool: /usr/bin/go test -benchmem -run=A$ \
-bench ABenchmarkGetPartsLength$ ieee-754/test

goos: linux
goarch: amd64
pkg: ieee-754/test
cpu: AMD Ryzen 5 5600H with Radeon Graphics

BenchmarkGetPartsLength-12   25748916   47.21 ns/op   0 B/op   0 allocs/op

PASS ok ieee-754/test 1.266s

Реализация алгоритмов

Теперь необходимо реализовать арифметические операции.

// Функция сложения.
func Sum(a, b float64) float64 {
	aFracLen := alt.GetPartsLength(a, 'f', -1, 64)
	bFracLen := alt.GetPartsLength(b, 'f', -1, 64)
	
	aModifier := math.Pow10(aFracLen)
	bModifier := math.Pow10(bFracLen)
	modifier := math.Max(aModifier, bModifier)
	
	a *= modifier
	b *= modifier
	
	return (a + b) / modifier
}

// Функция вычитания.
func Sub(a, b float64) float64 {
	aFracLen := alt.GetPartsLength(a, 'f', -1, 64)
	bFracLen := alt.GetPartsLength(b, 'f', -1, 64)
	
	aModifier := math.Pow10(aFracLen)
	bModifier := math.Pow10(bFracLen)
	
	modifier := math.Max(aModifier, bModifier)
	a *= modifier
	b *= modifier
	
	return (a - b) / modifier
}

// Функция умножения.
func Mul(a, b float64) float64 {
	aFracLen := alt.GetPartsLength(a, 'f', -1, 64)
	bFracLen := alt.GetPartsLength(b, 'f', -1, 64)
	
	aModifier := math.Pow10(aFracLen)
	bModifier := math.Pow10(bFracLen)
	
	modifier := aModifier * bModifier
	a *= aModifier
	b *= bModifier
	
	return a * b / modifier
}

// Функция деления.
func Div(a, b float64) float64 {
	aFracLen := alt.GetPartsLength(a, 'f', -1, 64)
	bFracLen := alt.GetPartsLength(b, 'f', -1, 64)

	aModifier := math.Pow10(aFracLen)
	bModifier := math.Pow10(bFracLen)

	minModifier := math.Min(aModifier, bModifier)
	maxModifier := math.Max(aModifier, bModifier)
	modifier := maxModifier / minModifier

	a *= aModifier
	b *= bModifier
	result := a / b

	if aModifier > bModifier {
		return result / modifier
	}

	if bModifier > aModifier {
		return result * modifier
	}

	return result
}

Наконец, давайте сравним по производительности обычную операцию с числами с плавающей запятой и безопасную:

func BenchmarkSafeMulFloat64(b *testing.B) {
	for i := 0; i < b.N; i++ {
		safecomp.Mul(2090.5, 8.61)
	}
}

func BenchmarkRegularMulFloat64(b *testing.B) {
	for i := 0; i < b.N; i++ {
		_ = 2090.5 * 8.61
	}
}

func BenchmarkMulDecimal(b *testing.B) {
	x := decimal.NewFromFloat(2090.5)
	у := decimal.NewFromFloat(8.61)

	for i := 0; i < b.N; i++ {
		x.Mul(y)
	}
}

Результат:

Running tool: /usr/bin/go test -benchmem -run=A$ \
-bench л(BenchmarkSafeMulFloat64|BenchmarkRegularMulFloat64| \
BenchmarkMulDecimal)$ github.com/zergon321/safecomp

goos: linux
goarch: amd64
pkg: github.com/zergon321/safecomp
cpu: AMD Ryzen 5 5600H with Radeon Graphics

BenchmarkSafeMulFloat64-12   16182450    74.24 ns/op   0 B/op   0 allocs/op
BenchmarkRegularMulFloat64-12 1000000000  0.2418 ns/op  0 B/op  0 allocs/op
BenchmarkMulDecimal-12      15515397    114.8 ns/op    80 B/op  2 allocs/op

PASS
ok github.com/zergon321/safecomp 3.408s

Это примерно в 313 раз медленнее обычной встроенной операции, но зато точно. Большая часть времени выполнения тратится на получение длины дробной части для обоих операндов. Возможно, это можно оптимизировать ещё больше. Также я поместил фрагмент кода умножения с самой популярной на Go реализацией Decimal. Как видите, безопасный метод вычисления float немного эффективнее (он дает те же результаты с //go:noinline).

Заключение

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