SBCL lead me on...

4 views
Skip to first unread message

Bob Felts

unread,
Sep 25, 2008, 8:52:59 PM9/25/08
to
Don Knuth, the author of "The Art of Computer Programming", has an
outstanding offer of a reward for the first report of an error in this
set of volumes. In exercise 4 of section 1.2.5, he states that log10
1000! = 2567.604_6_4.

SBCL (1.0.19.26, Mac OS X) says that

(log (loop for n from 1 to 1000 for f = 1 then (* f n)
finally (return f))
10d0)

is 2567.604_7_482272297

Could it be that Knuth was wrong?

Unfortunately, Maxima says:

(%i1) bfloat(log(1000!)/log(10));
(%o1) 2.567604_6_44222133b3


So Maxima agrees with Knuth and not SBCL.

Oh, well. So much for a chance at geek fame.

smallpond

unread,
Sep 26, 2008, 12:13:27 AM9/26/08
to


(loop for n from 1 to 1000 for f = 0d0 then (+ f (log n 10d0)) finally
(return f))
2567.6046442221304d0

--S

John Thingstad

unread,
Sep 26, 2008, 3:23:59 AM9/26/08
to

I think you should abide by the spirit of this claim. Knuth is assuming a
conventional language like C or Pascal (ALGOL?) with limited precision. If
you do the calculation here you will probably end up with 0. (The moral is
idiots get nothing.) Any how the question could be rephrased as what are
the number of digits in 1000! which is the average number of digits times
the number of digits. roughly 2.5 * 1000 = 2500. This is of course in
cents so about 2.5 dollars.

--------------
John Thingstad

Michael Weber

unread,
Sep 26, 2008, 5:16:31 AM9/26/08
to
On Sep 26, 6:13 am, smallpond <smallp...@juno.com> wrote:
>
> (loop for n from 1 to 1000 for f = 0d0 then (+ f (log n 10d0)) finally
> (return f))
> 2567.6046442221304d0

This is shorter:
(loop for n from 1 to 1000 summing (log n 10d0))

And here's one using SERIES:
(let ((range (scan-range :from 1 :upto 1000)))
(collect-sum (#Mlog range (series 10d0))))

--
Michael

Raymond Toy

unread,
Sep 26, 2008, 9:39:07 AM9/26/08
to
>>>>> "Bob" == Bob Felts <wr...@stablecross.com> writes:

Bob> Don Knuth, the author of "The Art of Computer Programming", has an
Bob> outstanding offer of a reward for the first report of an error in this
Bob> set of volumes. In exercise 4 of section 1.2.5, he states that log10
Bob> 1000! = 2567.604_6_4.

Bob> SBCL (1.0.19.26, Mac OS X) says that

Bob> (log (loop for n from 1 to 1000 for f = 1 then (* f n)
Bob> finally (return f))
Bob> 10d0)

Bob> is 2567.604_7_482272297

cmucl returns 2567.6046442221327d0.

Consider:

(/ (log (loop for n from 1 to 1000 for f = 1 then (* f n)
finally (return f)))
(log 10d0)) ->
2567.6047482272297d0

So, my guess is that sbcl is computing (log x 10d0) as (/ (log x) (log
10d0)). Since x is an integer, (log x) is a single-float. Hence,
only single-float accuracy for the result.

Ray

Waldek Hebisch

unread,
Sep 26, 2008, 11:23:34 AM9/26/08
to

Results form some other implementations:

sbcl 2567.6047482272297d0 (wrong)
Closure CL 2567.6047482272297D0 (wrong)
gcl Error: Can't print a non-number.
ecl Overflow
clisp Overflow
Poplog clisp Overflow

so, it seems that only cmucl gets it right.

> Oh, well. So much for a chance at geek fame.

ATM hitting bug in six different Lisp implenetations must do...

--
Waldek Hebisch
heb...@math.uni.wroc.pl

Kaz Kylheku

unread,
Sep 26, 2008, 5:02:57 PM9/26/08
to
On 2008-09-26, Bob Felts <wr...@stablecross.com> wrote:
> Don Knuth, the author of "The Art of Computer Programming", has an
> outstanding offer of a reward for the first report of an error in this
> set of volumes. In exercise 4 of section 1.2.5, he states that log10
> 1000! = 2567.604_6_4.
>
> SBCL (1.0.19.26, Mac OS X) says that
>
> (log (loop for n from 1 to 1000 for f = 1 then (* f n)
> finally (return f))
> 10d0)
>
> is 2567.604_7_482272297

In CLISP we get a floating point overflow: it doesn't like the large bignum
being coerced to a float.

But we can ``cheat'' using the equality

log (x/y) = log x - log y

i.e.
log x = log (x/y) + log y

Our bignum factorial has 2568 digits in it, so we can take y to be (expt 10
2568), and log y to be 2568.

We now end up taking the log of a bignum rational x/y, which CLISP does not
mind coercing to a float.

(+ (log (/ (loop for n from 1 to 1000

for f = 1 then (* f n) finally (return f))

(expt 10 2568))
10d0)
2568)

-> 2567.6046442221327d0

We can get more fractional digits by adding only 1 instead of 2568:

(+ (log (/ (loop for n from 1 to 1000

for f = 1 then (* f n) finally (return f))

(expt 10 2568))
10d0)
1)

-> 0.6046442221328487d0

Kaz Kylheku

unread,
Sep 26, 2008, 5:28:46 PM9/26/08
to
On 2008-09-26, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> Results form some other implementations:
>
> sbcl 2567.6047482272297d0 (wrong)
> Closure CL 2567.6047482272297D0 (wrong)
> gcl Error: Can't print a non-number.
> ecl Overflow
> clisp Overflow
> Poplog clisp Overflow
>
> so, it seems that only cmucl gets it right.
>
>> Oh, well. So much for a chance at geek fame.
>
> ATM hitting bug in six different Lisp implenetations must do...

The overflows are not bugs.

There is no requirement in the language spec that the log function must be able
to work with bignum integers which are outside of the range of the widest
floating point type.

The most desireable handling of the situation is to produce the right answer.
The second most desireable response is to signal an error. Producing the wrong
answer competes for third place with crashing.

Bob Felts

unread,
Sep 26, 2008, 7:44:16 PM9/26/08
to
Kaz Kylheku <kkyl...@gmail.com> wrote:

> On 2008-09-26, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> > Results form some other implementations:
> >
> > sbcl 2567.6047482272297d0 (wrong)
> > Closure CL 2567.6047482272297D0 (wrong)
> > gcl Error: Can't print a non-number.
> > ecl Overflow
> > clisp Overflow
> > Poplog clisp Overflow
> >
> > so, it seems that only cmucl gets it right.
> >
> >> Oh, well. So much for a chance at geek fame.
> >
> > ATM hitting bug in six different Lisp implenetations must do...
>
> The overflows are not bugs.
>
> There is no requirement in the language spec that the log function must be
> able to work with bignum integers which are outside of the range of the
> widest floating point type.

But why would anyone expect this? Surely I could use a Taylor series
even with bignums to calculate logs.

>
> The most desireable handling of the situation is to produce the right
> answer. The second most desireable response is to signal an error.
> Producing the wrong answer competes for third place with crashing.

It's interesting sbcl and Closure CL get the same wrong answer. Do they
have the same numerics package?

Rob Warnock

unread,
Sep 26, 2008, 9:14:45 PM9/26/08
to
Raymond Toy <raymo...@ericsson.com> wrote:
+---------------

| >>>>> "Bob" == Bob Felts <wr...@stablecross.com> writes:
| Bob> SBCL (1.0.19.26, Mac OS X) says that
| Bob> (log (loop for n from 1 to 1000 for f = 1 then (* f n)
| Bob> finally (return f))
| Bob> 10d0)
| Bob> is 2567.604_7_482272297
|
| cmucl returns 2567.6046442221327d0.
+---------------

Is that something that changed after CMUCL-19c?

cmu> (lisp-implementation-version)

"19c Release (19C)"
cmu> (log (loop for n from 1 to 1000 for f = 1 then (* f n)
finally (return f))
10d0)

2567.6047482272297d0
cmu>

Note that 19c does get an answer much closer to the one(s)
called "correct" if one sums the logs incrementally:

cmu> (loop for n from 1 to 1000 summing (log n 10d0))

2567.6046488768784d0
cmu>

but it's still different from yours. Finally, explicitly coercing
before the LOG gets still a third result, *almost* the same as
previously-posted "correct" answers:

cmu> (loop for n from 1 to 1000
summing (log (coerce n 'double-float) 10d0))

2567.60464422213d0
cmu>

???


-Rob

-----
Rob Warnock <rp...@rpw3.org>
627 26th Avenue <URL:http://rpw3.org/>
San Mateo, CA 94403 (650)572-2607

fvmar...@yahoo.com

unread,
Sep 27, 2008, 6:51:29 AM9/27/08
to
I don't have a copy of Knuth's "The Art of Computer Programming",
so I don't know what methods he discusses, but the standard numerical
approach to this problem would be to use Stirling's Approximation:

http://en.wikipedia.org/wiki/Stirling's_approximation

Here's a quick and dirty version, using the first four terms:

(defun stirling-approx (n)
(+ (* n (log n))
(- n)
(* 0.5d0 (log (* 2d0 pi n)))
(/ 1d0 12d0 n)))

(defun log10n! (n)
(/ (stirling-approx n) (log 10d0)))

(defun log10n!-error (n)
(/ (- 1d0) 360d0 (expt n 3) (log 10d0)))

CL-USER> (log10n! 1d3)
2567.6046442221336d0

CL-USER> (log10n!-error 1d3)
-1.206373560842366d-12

The absolute error in the formula in log10n! is at most 1.21e-12 for n
= 1000.

That said, I have no idea if log10n! is actually producing 16 good
digits.

-Frank

Scott Burson

unread,
Sep 27, 2008, 12:56:24 PM9/27/08
to
Aaargh!! I see this everywhere these days -- even in newspaper
stories! -- and while I'm sure that pointing it out here on c.l.l will
have no useful effect whatsoever except letting me blow off steam, I
can't resist.

The past tense of the verb "to lead" (long vowel) is spelled "led".
Yes, I know it's confusing, because it rhymes with the name of the
soft gray metal, which is spelled "lead" (short vowel). What's worse,
the past tense of "to read" (long vowel) is spelled "read" (short
vowel), because 'red' is another word altogether.

But you can deal with it anyway! Lots of us do! Today I lead,
yesterday you led, she has led many times in the past!

-- Scott

Raffael Cavallaro

unread,
Sep 27, 2008, 3:31:30 PM9/27/08
to
On 2008-09-27 12:56:24 -0400, Scott Burson <FSet...@gmail.com> said:

> The past tense of the verb "to lead" (long vowel) is spelled "led".
> Yes, I know it's confusing, because it rhymes with the name of the
> soft gray metal, which is spelled "lead" (short vowel). What's worse,
> the past tense of "to read" (long vowel) is spelled "read" (short
> vowel), because 'red' is another word altogether.

And here I thought the OP's subject line was an imperative exhortation
along the lines of:

"I for one welcome our new open source common lisp implementation
overlords; sbcl, lead me on!"


;^)

P.S. where does all this leave Led Zeppelin?

Bob Felts

unread,
Sep 27, 2008, 4:16:03 PM9/27/08
to
Scott Burson <FSet...@gmail.com> wrote:

You're absolutely right. I knew something didn't look right, but I
didn't follow through. My brain has been ruined by too much exposure to
lead (Pb) and too many hardware devices with LEDs. I repent in dust and
ashes and I appreciate you taking the time to hit me upside the head.

Pillsy

unread,
Sep 27, 2008, 4:35:03 PM9/27/08
to
On Sep 26, 5:28 pm, Kaz Kylheku <kkylh...@gmail.com> wrote:

> The most desireable handling of the situation is to produce the right answer.
> The second most desireable response is to signal an error.  Producing the wrong
> answer competes for third place with crashing.

The culprit in SBCL, as of version 10.1.17, appears to be the
following part of its LOG implementation:

(cond
((zerop base) 0f0) ; FIXME: type
((and (typep number '(integer (0) *))
(typep base '(integer (0) *)))
(coerce (/ (log2 number) (log2 base)) 'single-float))
(t (/ (log number) (log base))))

Taking the LOG of a bignum produces a single float, taking the LOG of
double produces a double float, and then the rules of type
contamination mean that the function returns the single float coerced
to a double float divided by a double float. The wrong double float.

Cheers,
Pillsy

Scott Burson

unread,
Sep 28, 2008, 12:22:29 AM9/28/08
to
On Sep 27, 1:16 pm, w...@stablecross.com (Bob Felts) wrote:

Thank you for taking my rant in a good spirit!

Cheers,
Scott

Nikodemus Siivola

unread,
Sep 28, 2008, 6:26:13 AM9/28/08
to
On Sep 26, 9:39 am, Raymond Toy <raymond....@ericsson.com> wrote:

> So, my guess is that sbcl is computing (logx 10d0) as (/ (logx) (log
> 10d0)).  Since x is an integer, (logx) is a single-float.  Hence,


> only single-float accuracy for the result.

Quite so. Fixed in 1.0.20.31.

Cheers,

-- Nikodemus

Scott Burson

unread,
Sep 28, 2008, 1:35:06 PM9/28/08
to
On Sep 27, 9:22 pm, Scott Burson <FSet....@gmail.com> wrote:
>
>
> Thank you for taking my rant in a good spirit!

... and don't be too hard on yourself. As I said, this error is
rampant these days; I'm sure you've picked it up from somewhere else.

-- Scott

Raymond Toy

unread,
Sep 29, 2008, 9:01:43 AM9/29/08
to
>>>>> "Rob" == Rob Warnock <rp...@rpw3.org> writes:

Rob> Raymond Toy <raymo...@ericsson.com> wrote:
Rob> +---------------
Rob> | >>>>> "Bob" == Bob Felts <wr...@stablecross.com> writes:
Rob> | Bob> SBCL (1.0.19.26, Mac OS X) says that
Rob> | Bob> (log (loop for n from 1 to 1000 for f = 1 then (* f n)
Rob> | Bob> finally (return f))
Rob> | Bob> 10d0)
Rob> | Bob> is 2567.604_7_482272297
Rob> |
Rob> | cmucl returns 2567.6046442221327d0.
Rob> +---------------

Rob> Is that something that changed after CMUCL-19c?

Yes. It seems to have been changed around somewhere around Jan,
2007. And mostly for these kinds of cases because I was
annoyed that (log <integer> 10d0) didn't have double-float accuracy.

I know the spec says (log number base) = (/ (log number) (log base))
which might imply only single precision results, but I didn't like
that. Basically, float contagion is applied before computing the
logs.

BTW, there's a typo in the spec. It says (log base number) = (/ (log
number) (log base)).

Rob> but it's still different from yours. Finally, explicitly coercing
Rob> before the LOG gets still a third result, *almost* the same as
Rob> previously-posted "correct" answers:

cmu> (loop for n from 1 to 1000

Rob> summing (log (coerce n 'double-float) 10d0))

Rob> 2567.60464422213d0

Because, previously (log n 10d0) was computed as (/ (log n) (log
10d0)), and (log n) is single-float, as required.

Ray

Thomas A. Russ

unread,
Oct 1, 2008, 12:31:13 PM10/1/08
to
Raffael Cavallaro <raffaelcavallaro@pas-d'espam-s'il-vous-plait-mac.com> writes:

> P.S. where does all this leave Led Zeppelin?

Wherever the aircraft that the zeppelin was following went?

--
Thomas A. Russ, USC/Information Sciences Institute

Raffael Cavallaro

unread,
Oct 1, 2008, 3:12:43 PM10/1/08
to
On 2008-10-01 12:31:13 -0400, t...@sevak.isi.edu (Thomas A. Russ) said:

> Wherever the aircraft that the zeppelin was following went?

meaning the Jefferson Airplane I suppose?

;^)

Kenny

unread,
Oct 1, 2008, 5:06:10 PM10/1/08
to

You mean Xah and White Rabbit, I guess.

kt

Reply all
Reply to author
Forward
0 new messages