saturate.nim 2.1 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980
  1. #
  2. #
  3. # The Nim Compiler
  4. # (c) Copyright 2012 Andreas Rumpf
  5. #
  6. # See the file "copying.txt", included in this
  7. # distribution, for details about the copyright.
  8. #
  9. ## Saturated arithmetic routines. XXX Make part of the stdlib?
  10. proc `|+|`*(a, b: BiggestInt): BiggestInt =
  11. ## saturated addition.
  12. result = a +% b
  13. if (result xor a) >= 0'i64 or (result xor b) >= 0'i64:
  14. return result
  15. if a < 0 or b < 0:
  16. result = low(result)
  17. else:
  18. result = high(result)
  19. proc `|-|`*(a, b: BiggestInt): BiggestInt =
  20. result = a -% b
  21. if (result xor a) >= 0'i64 or (result xor not b) >= 0'i64:
  22. return result
  23. if b > 0:
  24. result = low(result)
  25. else:
  26. result = high(result)
  27. proc `|abs|`*(a: BiggestInt): BiggestInt =
  28. if a != low(a):
  29. if a >= 0: result = a
  30. else: result = -a
  31. else:
  32. result = low(a)
  33. proc `|div|`*(a, b: BiggestInt): BiggestInt =
  34. # (0..5) div (0..4) == (0..5) div (1..4) == (0 div 4)..(5 div 1)
  35. if b == 0'i64:
  36. # make the same as ``div 1``:
  37. result = a
  38. elif a == low(a) and b == -1'i64:
  39. result = high(result)
  40. else:
  41. result = a div b
  42. proc `|mod|`*(a, b: BiggestInt): BiggestInt =
  43. if b == 0'i64:
  44. result = a
  45. else:
  46. result = a mod b
  47. proc `|*|`*(a, b: BiggestInt): BiggestInt =
  48. var
  49. resAsFloat, floatProd: float64
  50. result = a *% b
  51. floatProd = toBiggestFloat(a) # conversion
  52. floatProd = floatProd * toBiggestFloat(b)
  53. resAsFloat = toBiggestFloat(result)
  54. # Fast path for normal case: small multiplicands, and no info
  55. # is lost in either method.
  56. if resAsFloat == floatProd: return result
  57. # Somebody somewhere lost info. Close enough, or way off? Note
  58. # that a != 0 and b != 0 (else resAsFloat == floatProd == 0).
  59. # The difference either is or isn't significant compared to the
  60. # true value (of which floatProd is a good approximation).
  61. # abs(diff)/abs(prod) <= 1/32 iff
  62. # 32 * abs(diff) <= abs(prod) -- 5 good bits is "close enough"
  63. if 32.0 * abs(resAsFloat - floatProd) <= abs(floatProd):
  64. return result
  65. if floatProd >= 0.0:
  66. result = high(result)
  67. else:
  68. result = low(result)