Add `exact-integer-log` (similar to R7RS `exact-integer-sqrt`)

[jpellegrini]

exact-integer-log is to log what exact-integer-sqrt is to sqrt. It is actually easy to implement, since

• $\log_b(n) = \log_2(n) / \log_2(b)$.
• (integer-length n) $= \lfloor \log_2(n)\rfloor + 1$

We just need to put everything together, noting that $\log_b(n) = \log_2(n)/\log_2(b)$. Actually, we need a small adjustment, since

\lfloor \log_r(s)\rfloor / \lfloor \log_y(x)\rfloor
 \geq
\lfloor \log_r(s) / \log_y(x)\rfloor

(That is, it's not always an equality) - but that's easy to do.

Then we'll have this:

stklos> (exact-integer-log (+ (expt 7 20000) 5) 49)
10000
5

[jpellegrini]

• The Yacas computer algebra sysyem has something similar, but it only returns the integer part of the log. The function there is called IntLog: https://yacas.readthedocs.io/en/latest/reference_manual/programming.html?highlight=IntLog#IntLog

But the Yacas implementation is very slow.

[jpellegrini]

• The Yacas computer algebra sysyem has something similar

And Pari, too:

? logint(5*7^20000, 49)
%1 = 10000

But Pari is very fast.

[jpellegrini]

The "walking" at the end goes 1 at a time. It is possible to speed up this process, at the cost of some clarity:

(define (exact-integer-log n b)
  (unless (and (exact-integer? n) (positive? n))
    (error "log exponent ~S is not a positive exact integer" n))
  (unless (and (exact-integer? b) (> b 1))
    (error "log base ~S is not exact integer greater than one" b))
  ;; This is easier and faster than R7RS exact-integer-sqrt, since:
  ;; * log(n b) = log_2(n) / log_2(b).
  ;; * integer-length(n) = floor(log_2(n)) + 1
  ;; We just put it all together.
  (let ((log2-n (- (integer-length n) 1))
        (log2-b (- (integer-length b) 1)))
    (let ((L (quotient log2-n log2-b))
          (step 1))
      ;; Ok, so FLOOR(a/b) can be smaller than FLOOR(a)/FLOOR(b).
      ;; That means we may have overestimated L -- but we should be
      ;; close. If we overestimated, then (expt b L) will be larger
      ;; than n, and we can subtract from it until we find the correct
      ;; number.
      ;; We do increase the step at each iteration in order to speed
      ;; up the process. If we go too far, we correct it later. This is
      ;; fast.
      (while (> (expt b L) n)
        (set! L (- L step))
        (set! step (+ step 2)))

      ;; If we overestimated, then we had to walk down from L.
      ;; In this case, since we added 2 to the step at each
      ;; iteration, we may now have UNDERestimated. But then,
      ;; we undo the last iteration (the "go back" line below),
      ;; and start walking again, this time with a fixed step
      ;; equal to 1.
      (when (> step 1)            ; did we need to walk?
        (set! L (- (+ L step) 2)) ; go back...
        (while (> (expt b L) n)
          (set! L (- L 1))))      ; don't grow the step this time!

      (values L (- n (expt b L))))))

I can add this to the PR if you think it's OK @egallesio

[jpellegrini]

I can add this to the PR if you think it's OK @egallesio

Added as a separate commit, to make things easier.

[egallesio]

Sorry for the long time to merge this one. Finally, it's done. Thanks.