Skip to content

Instantly share code, notes, and snippets.

@Bike
Last active May 8, 2023 00:31
Show Gist options
  • Save Bike/a53d13430aa2adcfeb05cf7d2be9724c to your computer and use it in GitHub Desktop.
Save Bike/a53d13430aa2adcfeb05cf7d2be9724c to your computer and use it in GitHub Desktop.
Notes on optimizing CL arithmetic

Bitwise operations on signed magnitude bignums

Because CL specifies bitwise operations treat integers as being in two's complement, if your bignums use a signed magnitude representation they are nontrivial. You will have to test signs while doing your arithmetic operations, and do different things depending. Generally, keep in mind that for integer x, (- x) = (lognot (1- x)) = (1+ (lognot x)), and that the number of bignum elements in the output may exceed those in the input. (Consider (logand -2 -3): It's -4, and 4 needs one more bit to represent than 2 or 3 do.)

For bitwise operations of two operands, you can understand the sign of the result by thinking of the infinite left bits. For example, if logior is given one positive and one negative operand, the result is negative, because for each of these infinitely many bits, (logior 0 1) = 1.

Here are some particular identities that are useful, derived from basic boolean algebra plus the above. (- (foo ...)) indicates that the result of (foo ...) will be negative, and so (- (foo ...)) is the absolute magnitude of this result. So for example, (- (logand (- x) (- y))) = (1+ (logior (1- x) (1- y))) indicates that if logand is given two negative arguments, the result is a negative number with magnitude (1+ (logior (1- x) (1- y))).

  • (integer-length (- x)) = (integer-length (1- x))
  • (logand x (- y)) = (logandc2 x (1- y))
  • (- (logand (- x) (- y))) = (1+ (logior (1- x) (1- y)))
  • (logandc2 x (- y)) = (logand x (1- y))
  • (- (logandc2 (- x) y)) = (1+ (logior (1- x) y))
  • (logandc2 (- x) (- y)) = (logandc1 (1- x) (1- y))
  • (- (logeqv x y)) = (1+ (logxor x y))
  • (logeqv x (- y)) = (logxor x (1- y))
  • (- (logeqv (- x) (- y))) = (1+ (logxor (1- x) (1- y)))
  • (- (logior x (- y))) = (1+ (logandc1 x (1- y)))
  • (- (logior (- x) (- y))) = (1+ (logand (1- x) (1- y)))
  • (- (lognand x y)) = (1+ (logand x y))
  • (- (lognand x (- y))) = (1+ (logandc2 x (1- y)))
  • (lognand (- x) (- y)) = (logior (1- x) (1- y))
  • (- (lognor x y)) = (1+ (logior x y))
  • (lognor x (- y)) = (logandc1 x (1- y))
  • (lognor (- x) (- y)) = (logand (1- x) (1- y))
  • (- (logorc2 x y)) = (1+ (logandc1 x y))
  • (logorc2 x (- y)) = (logior x (1- y))
  • (- (logorc2 (- x) y)) = (1+ (logand (1- x) y))
  • (- (logorc2 (- x) (- y))) = (1+ (logandc2 (1- x) (1- y)))
  • (- (logxor x (- y))) = (1+ (logxor x (1- y)))
  • (logxor (- x) (- y)) = (logxor (1- x) (1- y))
  • (logcount (- x)) = (logcount (lognot (1- x))) = (integer-length (1- x))

ash also needs some special casing if the integer being shifted is a negative bignum being shifted to the right. For a negative, (- (ash (- x) a)) = (1+ (ash (1- x) a)).

Comparisons

For the following, cmp is a comparison operator, i.e. < or whatever.

Integer-to-integer comparisons

Since bignums are out of the range of fixnums, (> bignum fixnum) = (plusp bignum) and so on. And when comparing bignums you can obviously start by comparing their length in words/whatever.

Integer-to-rational comparisons

Obviously (cmp a x/y) = (cmp (* y a) x), but (* y a) may be a bignum requiring some consing. You may avoid this by using these identities (cribbed from SBCL): (< a x/y) = (< a (ceiling x y)), (< x/y z) = (< (floor x y) a), (> a x/y) = (> a (floor x y)), (> x/y a) = (> (ceiling x y) a).

Float-to-rational comparisons

CLHS 12.1.4.1 specifies that if f is a float and r a rational, (cmp f r) = (cmp (rational f) r). But this might mean consing up a ratio, which is bad.

For an integer, we can instead check the integer part of the float against the integer. (< (truncate f) i) implies (< f i) and (> (truncate f) i) implies (> f i) regardless of signs. When (= (truncate f) i), you have to check the fractional part of f. There are at least two ways to do this. If (> f (truncate f)) or (> f 0), (> f i), and if (< f (truncate f)) or (< f 0), (< f i), and otherwise (= f i).

Not sure about ratios yet.

Rationalization

cl:rational can be done very quickly if your float radix is 2 and you're on a binary machine (and if you're not, you need more help than I can give you). Say the magnitude of your float f is i * 2^e with i, e integral (what you get out of integer-decode-float). In general, (abs (rational f)) is then (* i (expt 2 e)), but you can skip most of the general computation. If e is positive, this is obviously just (ash i e). If e is negative, find the largest g such that i/2^g is still an integer, or put another way, find the least significant 1 bit of i (most processors have an instruction for this, "find first set"). If -e < g, your result is again (ash i e) and you're only shifting out zero bits. If not, the result is a ratio i/2^-e. But since you know the denominator is a power of two, you don't need to compute a gcd or anything to reduce the ratio: just return the ratio with numerator (ash i (- g)) and denominator (ash 1 (- (- e) g)), and you know it's in lowest terms.

In code,

(defun %rational (signif exponent)
  (if (minusp exponent)
      (let ((pexp (- exponent))
            (g (find-first-set signif)))
        (if (> pexp g)
            (%make-ratio (ash signif (- g)) (ash 1 (- pexp set)))
            (ash signif exponent)))
      (ash signif exponent)))

(defun rational (f)
  (etypecase f
    (rational f)
    (float
     (multiple-value-bind (significand exponent sign)
         (integer-decode-float f)
       (if (= significand 0)
           0
           (let ((abs (%rational significand exponent)))
             (if (= sign -1)
                 (- abs)
                 abs)))))))
@moon-chilled
Copy link

(< float rational) can be done with a plain ceiling division, after checking exponent and scaling numerator/denominator so the result is 53 bits (or 24 for singles). In particular, note that float can only be = to rational if the denominator of the latter is a power of two, which is easy to check for and rule out before anything else.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment