First power of 2 that has leading decimal digits of 12 in Common Lisp

Question

In the webpage of Rosetta Code (http://rosettacode.org/wiki/First_power_of_2_that_has_leading_decimal_digits_of_12#C) there is the following challenge: Find the first power of 2 that has leading decimal digits of 12.

I have written several versions in Common Lisp, but I fail miserably in terms of performance compared to what other languages report.

The following code is one of the many versions that I have tried.

(defun log10 (x) (/ (log x) (log 10)))

(defun first-digits (n l)
    (let* ((len-n (1+ (floor (log10 n))))
           (tens (expt 10 (- len-n l))) )
       (truncate n tens) ))


(defun p (rem n)
    (do* ((len-rem (1+ (floor (log10 rem))))
                (i 0 (1+ i)) 
        (k 1 (* 2 k)) )
       ((= n 0) (1- i))
    (when (= rem (first-digits k len-rem))
        (decf n) )))

The performance is really poor, but I refuse to admit that Common Lisp is slower than any competitors. Any idea of how to achieve the run time of a few seconds reported by C#, C++, etc.?

Answer 1

Here is an answer which is now substantially equivalent to yours (based on the idea in this maths SE answer, but a little more general and a little more optimized. In particular:

• it lets you specify both the number you want to raise to a power;
• it lets you specify the base you want to work in (so not just 10);
• it knows that (log x b) is the log of x in base b;
• it computes the log you need just once;
• it reorganizes some of the float operations to avoid overflows (at some point your version will call (expt 10 n) with a big enough n thay you'll get a floating point overflow.

I think (in fact I'm reasonably sure) that the remaining thing that is making it slow is float consing. It ought to be possible to prevail upon a compiler to avoid this, and I spent a small amount of time trying but to no avail.

Anyway, here it is.

(defun p (L n &key (b 2) (base 10.0d0))
  ;; nth occurence of power of b with leading digits L in base:
  ;; returns power (p below)
  (let ((k-1 (floor (log L base)))      ;(digits of L in base) - 1
        (logb (log (float b 1.0d0) base))) ;log of b in base
    (do* ((p 0 (1+ p))                  ;power
          (d 1 (floor (expt base (+ (rem (* p logb) 1) k-1)))) ;digits
          (c (if (= L 1) 1 0)            ;count of hits
             (if (= L d) (1+ c) c)))
         ((= c n) p))))

---

Below is an earlier, useless answer: I've left it for posterity.

Note that the first part what's below deals with only the two-leading-digit case because I didn't bother looking at the original rosetta code thing: see the end for a generalisation.

First of all to pull out the first two decimal digits of something, you want to start by divide it by a suitable power of 10 that the result is a two digit number. That means taking logs, but you can cheat: if you know the number of bits, b, in the number then the number is 2^b + change (assuming it's a positive integer), and knowing the number of bits is very very quick. And then if 10^x = 2^y then x = y/(log_2 10) where log_2 is log base 2. And this is a constant.

So you can write this function after fiddling around with pen and paper for a bit:

(defun leading-two-digits (n)
  (truncate (truncate n (expt 10 (1- (truncate (integer-length n)
                                               (load-time-value (log 10 2))))))
            10))

Note the use of load-time-value to compute the (log 10 2) just once. And note also that this relies on either having correct integer arithmetic or at least a function which will tell you what the number of bits in a number would be if you did have.

So now

> (leading-two-decimal-digits 59468049869823987435)
5
9

OK, looks good (extensive testing...).

So now you just need to start from 1 and successively multiply by the base, looking for the leading two digits being what you care about. As a hack if the base is 2 you can just left shift: my original function assumed the base was 2 and always left shifted, this one has a configurable base but still special-cases 2, which I suspect may no longer help at all.

(defun nth-ld-power (n &key (b 2) (leading 1) (second 2))
  ;; iterate is tfeb.github.io/tfeb-lisp-hax/, also should be in
  ;; Quicklisp
  (iterate looking ((v 1)               ;value
                    (p 0)               ;power (so we can report it)
                    (c 0))              ;hit count
    (multiple-value-bind (d1 d2) (leading-two-decimal-digits v)
      (if (and (= d1 leading) (= d2 second))
          (let ((c+ (1+ c)))
            (if (= c+ n)
                (values v p)
              (looking (if (= b 2) (ash v 1) (* v b)) (1+ p) c+)))
        (looking (if (= b 2) (ash v 1) (* v b)) (1+ p) c)))))

Again you would want to deal with the multiple-leading-digits case by passing a list of leading digits, but.

OK, so now:

> (time (nth-ld-power 2))
Timing the evaluation of (nth-ld-power 2)

User time    =        0.000
System time  =        0.000
Elapsed time =        0.000
Allocation   = 2848 bytes
0 Page faults
GC time      =        0.000
1208925819614629174706176
80

> (time (nth-ld-power 10))
Timing the evaluation of (nth-ld-power 10)

User time    =        0.000
System time  =        0.000
Elapsed time =        0.000
Allocation   = 40704 bytes
0 Page faults
GC time      =        0.000
124330809102446660538845562036705210025114037699336929360115994223289874253133343883264
286

OK, this is too quick to measure, let's make it do some real work:

> (time (nth-ld-power 1000))
Timing the evaluation of (nth-ld-power 1000)

User time    =        3.264
System time  =        0.039
Elapsed time =        3.384
Allocation   = 845193424 bytes
517 Page faults
GC time      =        0.020
12[...long number truncated here ...]
28745

So that's at least reasonable performance, I think. These figures were for LispWorks on a 2013 macbook. Obviously you can rewrite it without Tim's iterate, but I like it.

---

Here is a generalisation to the n-digit case.

(defun leading-decimal-digits (n m)
  ;; Return M leading digits from N: no sanity check
  (iterate peel ((q (truncate n (expt 10 (- (truncate (integer-length n)
                                                      (load-time-value (log 10 2)))
                                            (- m 1)))))
                 (digits '()))
    (multiple-value-bind (qq d) (truncate q 10)
      (cond
       ((zerop qq)
        (cons d digits))
       ((< qq 10)
         (list* qq d digits))
       (t
        (peel qq (cons d digits)))))))

(defun nth-power-with-leading-decimal-digits (n b leading)
  ;; Find the N'th occurrence of a power of B whose leading digits are
  ;; LEADING (a list of digits).  Return the value of B^P and P, the
  ;; power.
  (let ((nleading (length leading)))
    (iterate looking ((v 1)               ;value
                      (p 0)               ;power (so we can report it)
                      (c 0))              ;hit count
      (if (equal (leading-decimal-digits v nleading) leading)
          (let ((c+ (1+ c)))
            (if (= c+ n)
                (values v p)
              (looking (if (= b 2) (ash v 1) (* v b)) (1+ p) c+)))
        (looking (if (= b 2) (ash v 1) (* v b)) (1+ p) c)))))

So:

> (nth-power-with-leading-decimal-digits 2 2 '(1 2 8))
12855504354071922204335696738729300820177623950262342682411008
203

> (nth-power-with-leading-decimal-digits 2 3 '(2 7 8))
278954761343915929031866324148580803686773879062609352173281933430969939572023921519256921927359964084535215107090906022143908601839272147120823008337941481521208646465304746378648054338849857759629806700446921838039313884792762356955010344065298744426691826196079923777796821513452648753573059469525738664313409324728161550430310432705576201607066435772343529511415605105251217669677767280155388600975280964482318641251059803074701681895639266095014303466041595626938522373004626313927779288797843485898785327074755040298312905780373591994824646107875196292150692772609024142420597241695173176699988995274347532223981851419958515807471788303361043192244700492703668505222371093498892685289816397935636213340656919509760640657215827785179171568543835684943671352357398679232652259639784388841436515656182938820127024910154405205830179436914341715546500034485031465189386918271898229029201
1860

2^203 is the second occurrence of a power of 2 beginning with 128, and 3^1860 is the second occurrence of a power of 3 beginning 278. The second one took 0.017 seconds.