1 /++ 2 Note: 3 The module doesn't provide full arithmetic API for now. 4 +/ 5 module mir.bignum.fp; 6 7 import std.traits; 8 import mir.bitop; 9 import mir.utility; 10 11 package enum half(size_t hs) = (){ 12 import mir.bignum.fixed: UInt; 13 UInt!hs ret; ret.signBit = true; return ret; 14 }(); 15 16 /++ 17 Software floating point number. 18 19 Params: 20 coefficientSize = coefficient size in bits 21 22 Note: the implementation doesn't support NaN and Infinity values. 23 +/ 24 struct Fp(size_t coefficientSize, Exp = sizediff_t) 25 if ((is(Exp == int) || is(Exp == long)) && coefficientSize % (size_t.sizeof * 8) == 0 && coefficientSize >= (size_t.sizeof * 8)) 26 { 27 import mir.bignum.fixed: UInt; 28 29 bool sign; 30 Exp exponent; 31 UInt!coefficientSize coefficient; 32 33 /++ 34 +/ 35 nothrow 36 this(bool sign, Exp exponent, UInt!coefficientSize normalizedCoefficient) 37 { 38 this.coefficient = normalizedCoefficient; 39 this.exponent = exponent; 40 this.sign = sign; 41 } 42 43 /++ 44 Constructs $(LREF Fp) from hardaware floating point number. 45 Params: 46 value = Hardware floating point number. Special values `nan` and `inf` aren't allowed. 47 normalize = flag to indicate if the normalization should be performed. 48 +/ 49 this(T)(const T value, bool normalize = true) 50 @safe pure nothrow @nogc 51 if (isFloatingPoint!T && T.mant_dig <= coefficientSize) 52 { 53 import mir.math.common : fabs; 54 import mir.math.ieee : frexp, signbit, ldexp; 55 assert(value == value); 56 assert(value.fabs < T.infinity); 57 this.sign = value.signbit != 0; 58 if (value == 0) 59 return; 60 T x = value.fabs; 61 int exp; 62 { 63 enum scale = T(2) ^^ T.mant_dig; 64 x = frexp(x, exp) * scale; 65 exp -= T.mant_dig; 66 } 67 static if (T.mant_dig < 64) 68 { 69 this.coefficient = UInt!coefficientSize(cast(ulong)cast(long)x); 70 } 71 else 72 static if (T.mant_dig == 64) 73 { 74 this.coefficient = UInt!coefficientSize(cast(ulong)x); 75 } 76 else 77 { 78 enum scale = T(2) ^^ 64; 79 enum scaleInv = 1 / scale; 80 x *= scaleInv; 81 long high = cast(long) x; 82 if (high > x) 83 --high; 84 x -= high; 85 x *= scale; 86 this.coefficient = (UInt!coefficientSize(ulong(high)) << 64) | cast(ulong)x; 87 } 88 if (normalize) 89 { 90 auto shift = cast(int)this.coefficient.ctlz; 91 exp -= shift; 92 this.coefficient <<= shift; 93 } 94 else 95 { 96 int shift = T.min_exp - T.mant_dig - exp; 97 if (shift > 0) 98 { 99 this.coefficient >>= shift; 100 exp = T.min_exp - T.mant_dig; 101 } 102 } 103 this.exponent = exp; 104 } 105 106 static if (coefficientSize == 128) 107 /// 108 version(mir_bignum_test) 109 @safe pure @nogc nothrow 110 unittest 111 { 112 enum h = -33.0 * 2.0 ^^ -10; 113 auto f = Fp!64(h); 114 assert(f.sign); 115 assert(f.exponent == -10 - (64 - 6)); 116 assert(f.coefficient == 33UL << (64 - 6)); 117 assert(cast(double) f == h); 118 119 // CTFE 120 static assert(cast(double) Fp!64(h) == h); 121 122 f = Fp!64(-0.0); 123 assert(f.sign); 124 assert(f.exponent == 0); 125 assert(f.coefficient == 0); 126 127 // subnormals 128 static assert(cast(float) Fp!64(float.min_normal / 2) == float.min_normal / 2); 129 static assert(cast(float) Fp!64(float.min_normal * float.epsilon) == float.min_normal * float.epsilon); 130 // subnormals 131 static assert(cast(double) Fp!64(double.min_normal / 2) == double.min_normal / 2); 132 static assert(cast(double) Fp!64(double.min_normal * double.epsilon) == double.min_normal * double.epsilon); 133 // subnormals 134 static assert(cast(real) Fp!64(real.min_normal / 2) == real.min_normal / 2); 135 static assert(cast(real) Fp!64(real.min_normal * real.epsilon) == real.min_normal * real.epsilon); 136 137 enum d = cast(float) Fp!64(float.min_normal / 2, false); 138 139 // subnormals 140 static assert(cast(float) Fp!64(float.min_normal / 2, false) == float.min_normal / 2, d.stringof); 141 static assert(cast(float) Fp!64(float.min_normal * float.epsilon, false) == float.min_normal * float.epsilon); 142 // subnormals 143 static assert(cast(double) Fp!64(double.min_normal / 2, false) == double.min_normal / 2); 144 static assert(cast(double) Fp!64(double.min_normal * double.epsilon, false) == double.min_normal * double.epsilon); 145 // subnormals 146 static assert(cast(real) Fp!64(real.min_normal / 2, false) == real.min_normal / 2); 147 static assert(cast(real) Fp!64(real.min_normal * real.epsilon, false) == real.min_normal * real.epsilon); 148 } 149 150 static if (coefficientSize == 128) 151 /// Without normalization 152 version(mir_bignum_test) 153 @safe pure @nogc nothrow 154 unittest 155 { 156 auto f = Fp!64(-33.0 * 2.0 ^^ -10, false); 157 assert(f.sign); 158 assert(f.exponent == -10 - (double.mant_dig - 6)); 159 assert(f.coefficient == 33UL << (double.mant_dig - 6)); 160 } 161 162 /++ 163 +/ 164 this(size_t size)(UInt!size integer, bool normalizedInteger = false) 165 nothrow 166 { 167 import mir.bignum.fixed: UInt; 168 static if (size < coefficientSize) 169 { 170 if (normalizedInteger) 171 { 172 this(false, Exp(size) - coefficientSize, integer.rightExtend!(coefficientSize - size)); 173 } 174 else 175 { 176 this(integer.toSize!coefficientSize, false); 177 } 178 } 179 else 180 { 181 this.exponent = size - coefficientSize; 182 if (!normalizedInteger) 183 { 184 auto c = integer.ctlz; 185 integer <<= c; 186 this.exponent -= c; 187 } 188 static if (size == coefficientSize) 189 { 190 coefficient = integer; 191 } 192 else 193 { 194 enum N = coefficient.data.length; 195 version (LittleEndian) 196 coefficient.data = integer.data[$ - N .. $]; 197 else 198 coefficient.data = integer.data[0 .. N]; 199 enum tailSize = size - coefficientSize; 200 auto cr = integer.toSize!tailSize.opCmp(half!tailSize); 201 if (cr > 0 || cr == 0 && coefficient.bt(0)) 202 { 203 if (auto overflow = coefficient += 1) 204 { 205 coefficient = half!coefficientSize; 206 exponent++; 207 } 208 } 209 } 210 } 211 } 212 213 static if (coefficientSize == 128) 214 /// 215 version(mir_bignum_test) 216 @safe pure @nogc 217 unittest 218 { 219 import mir.bignum.fixed: UInt; 220 221 auto fp = Fp!128(UInt!128.fromHexString("afbbfae3cd0aff2714a1de7022b0029d")); 222 assert(fp.exponent == 0); 223 assert(fp.coefficient == UInt!128.fromHexString("afbbfae3cd0aff2714a1de7022b0029d")); 224 225 fp = Fp!128(UInt!128.fromHexString("afbbfae3cd0aff2714a1de7022b0029d"), true); 226 assert(fp.exponent == 0); 227 assert(fp.coefficient == UInt!128.fromHexString("afbbfae3cd0aff2714a1de7022b0029d")); 228 229 fp = Fp!128(UInt!128.fromHexString("ae3cd0aff2714a1de7022b0029d")); 230 assert(fp.exponent == -20); 231 assert(fp.coefficient == UInt!128.fromHexString("ae3cd0aff2714a1de7022b0029d00000")); 232 233 fp = Fp!128(UInt!128.fromHexString("e7022b0029d")); 234 assert(fp.exponent == -84); 235 assert(fp.coefficient == UInt!128.fromHexString("e7022b0029d000000000000000000000")); 236 237 fp = Fp!128(UInt!64.fromHexString("e7022b0029d")); 238 assert(fp.exponent == -84); 239 assert(fp.coefficient == UInt!128.fromHexString("e7022b0029d000000000000000000000")); 240 241 fp = Fp!128(UInt!64.fromHexString("e7022b0029dd0aff"), true); 242 assert(fp.exponent == -64); 243 assert(fp.coefficient == UInt!128.fromHexString("e7022b0029dd0aff0000000000000000")); 244 245 fp = Fp!128(UInt!64.fromHexString("e7022b0029d")); 246 assert(fp.exponent == -84); 247 assert(fp.coefficient == UInt!128.fromHexString("e7022b0029d000000000000000000000")); 248 249 fp = Fp!128(UInt!192.fromHexString("ffffffffffffffffffffffffffffffff1000000000000000")); 250 assert(fp.exponent == 64); 251 assert(fp.coefficient == UInt!128.fromHexString("ffffffffffffffffffffffffffffffff")); 252 253 fp = Fp!128(UInt!192.fromHexString("ffffffffffffffffffffffffffffffff8000000000000000")); 254 assert(fp.exponent == 65); 255 assert(fp.coefficient == UInt!128.fromHexString("80000000000000000000000000000000")); 256 257 fp = Fp!128(UInt!192.fromHexString("fffffffffffffffffffffffffffffffe8000000000000000")); 258 assert(fp.exponent == 64); 259 assert(fp.coefficient == UInt!128.fromHexString("fffffffffffffffffffffffffffffffe")); 260 261 fp = Fp!128(UInt!192.fromHexString("fffffffffffffffffffffffffffffffe8000000000000001")); 262 assert(fp.exponent == 64); 263 assert(fp.coefficient == UInt!128.fromHexString("ffffffffffffffffffffffffffffffff")); 264 } 265 266 /// 267 ref Fp opOpAssign(string op : "*")(Fp rhs) nothrow return 268 { 269 this = this.opBinary!op(rhs); 270 return this; 271 } 272 273 /// 274 Fp opBinary(string op : "*")(Fp rhs) nothrow const 275 { 276 return cast(Fp) .extendedMul(this, rhs); 277 } 278 279 static if (coefficientSize == 128) 280 /// 281 version(mir_bignum_test) 282 @safe pure @nogc 283 unittest 284 { 285 import mir.bignum.fixed: UInt; 286 287 auto a = Fp!128(0, -13, UInt!128.fromHexString("dfbbfae3cd0aff2714a1de7022b0029d")); 288 auto b = Fp!128(1, 100, UInt!128.fromHexString("e3251bacb112c88b71ad3f85a970a314")); 289 auto fp = a * b; 290 assert(fp.sign); 291 assert(fp.exponent == 100 - 13 + 128); 292 assert(fp.coefficient == UInt!128.fromHexString("c6841dd302415d785373ab6d93712988")); 293 } 294 295 /// 296 T opCast(T)() nothrow const 297 if (is(Unqual!T == bool)) 298 { 299 return coefficient != 0; 300 } 301 302 /// 303 T opCast(T, bool noHalf = false)() nothrow const 304 if (isFloatingPoint!T) 305 { 306 import mir.math.ieee: ldexp; 307 auto exp = cast()exponent; 308 static if (coefficientSize == 32) 309 { 310 Unqual!T c = cast(uint) coefficient; 311 } 312 else 313 static if (coefficientSize == 64) 314 { 315 Unqual!T c = cast(ulong) coefficient; 316 } 317 else 318 { 319 enum shift = coefficientSize - T.mant_dig; 320 enum rMask = (UInt!coefficientSize(1) << shift) - UInt!coefficientSize(1); 321 enum rHalf = UInt!coefficientSize(1) << (shift - 1); 322 enum rInc = UInt!coefficientSize(1) << shift; 323 UInt!coefficientSize adC = coefficient; 324 static if (!noHalf) 325 { 326 auto cr = (coefficient & rMask).opCmp(rHalf); 327 if ((cr > 0) | (cr == 0) & coefficient.bt(shift)) 328 { 329 if (auto overflow = adC += rInc) 330 { 331 adC = half!coefficientSize; 332 exp++; 333 } 334 } 335 } 336 adC >>= shift; 337 exp += shift; 338 Unqual!T c = cast(ulong) adC; 339 static if (T.mant_dig > 64) // 340 { 341 static assert (T.mant_dig <= 128); 342 c += ldexp(cast(T) cast(ulong) (adC >> 64), 64); 343 } 344 } 345 if (sign) 346 c = -c; 347 static if (exp.sizeof > int.sizeof) 348 { 349 import mir.utility: min, max; 350 exp = exp.max(int.min).min(int.max); 351 } 352 return ldexp(c, cast(int)exp); 353 } 354 355 static if (coefficientSize == 128) 356 /// 357 version(mir_bignum_test) 358 @safe pure @nogc 359 unittest 360 { 361 import mir.bignum.fixed: UInt; 362 auto fp = Fp!128(1, 100, UInt!128.fromHexString("e3251bacb112cb8b71ad3f85a970a314")); 363 assert(cast(double)fp == -0xE3251BACB112C8p+172); 364 } 365 366 static if (coefficientSize == 128) 367 /// 368 version(mir_bignum_test) 369 @safe pure @nogc 370 unittest 371 { 372 import mir.bignum.fixed: UInt; 373 auto fp = Fp!128(1, 100, UInt!128.fromHexString("e3251bacb112cb8b71ad3f85a970a314")); 374 static if (real.mant_dig == 64) 375 assert(cast(real)fp == -0xe3251bacb112cb8bp+164L); 376 } 377 378 static if (coefficientSize == 128) 379 /// 380 version(mir_bignum_test) 381 @safe pure @nogc 382 unittest 383 { 384 import mir.bignum.fixed: UInt; 385 auto fp = Fp!64(1, 100, UInt!64(0xe3251bacb112cb8b)); 386 version (DigitalMars) 387 { 388 // https://issues.dlang.org/show_bug.cgi?id=20963 389 assert(cast(double)fp == -0xE3251BACB112C8p+108 390 || cast(double)fp == -0xE3251BACB112D0p+108); 391 } 392 else 393 { 394 assert(cast(double)fp == -0xE3251BACB112C8p+108); 395 } 396 } 397 // -0x1.c64a375962259p+163 = 398 // -0xe.3251bacb112cb8bp+160 = 399 // -0x1.c64a37596225ap+163 = 400 // -0xe.3251bacb112cb8bp+160 = 401 static if (coefficientSize == 128) 402 /// 403 version(mir_bignum_test) 404 @safe pure @nogc 405 unittest 406 { 407 import mir.bignum.fixed: UInt; 408 auto fp = Fp!64(1, 100, UInt!64(0xe3251bacb112cb8b)); 409 static if (real.mant_dig == 64) 410 assert(cast(real)fp == -0xe3251bacb112cb8bp+100L); 411 } 412 413 /// 414 T opCast(T : Fp!newCoefficientSize, size_t newCoefficientSize)() nothrow const 415 { 416 auto ret = Fp!newCoefficientSize(coefficient, true); 417 ret.exponent += exponent; 418 ret.sign = sign; 419 return ret; 420 } 421 422 static if (coefficientSize == 128) 423 /// 424 version(mir_bignum_test) 425 @safe pure @nogc 426 unittest 427 { 428 import mir.bignum.fixed: UInt; 429 auto fp = cast(Fp!64) Fp!128(UInt!128.fromHexString("afbbfae3cd0aff2784a1de7022b0029d")); 430 assert(fp.exponent == 64); 431 assert(fp.coefficient == UInt!64.fromHexString("afbbfae3cd0aff28")); 432 } 433 } 434 435 /// 436 Fp!(coefficientizeA + coefficientizeB) extendedMul(size_t coefficientizeA, size_t coefficientizeB)(Fp!coefficientizeA a, Fp!coefficientizeB b) 437 @safe pure nothrow @nogc 438 { 439 import mir.bignum.fixed: extendedMul; 440 auto coefficient = extendedMul(a.coefficient, b.coefficient); 441 auto exponent = a.exponent + b.exponent; 442 auto sign = a.sign ^ b.sign; 443 if (!coefficient.signBit) 444 { 445 --exponent; 446 coefficient = coefficient.smallLeftShift(1); 447 } 448 return typeof(return)(sign, exponent, coefficient); 449 }