Prev: Cubic equations and companion matrices
Next: Why does JH Conway say -1 is prime?
From: Richard Tobin on 6 Feb 2010 12:25
In article <20100206120206.850$lA(a)newsreader.com>,
David W. Cantrell <DWCantrell(a)sigmaxi.net> wrote:

>> Could you determine whether it fails for 3121515!. I don't have the
>> facilities to hand to check it accurately.

That was just the first case where 64-bit IEEE floating point numbers
give different results, calculating the factorial by adding the logs.
Of course, I couldn't tell whether it was the accumulated floating
point error or a real failure.

>Any other suggestions?

Using "long double", the first candidates are 48655817 and 91948153
(beyond that they occur often, suggesting that the floating point
error has become significant).

-- Richard
--
Please remember to mention me / in tapes you leave behind.
From: David W. Cantrell on 6 Feb 2010 12:33
Chip Eastham <hardmath(a)gmail.com> wrote:
> On Feb 6, 1:56=A0am, David W. Cantrell <DWCantr...(a)sigmaxi.net> wrote:
> > "David W. Cantrell" <DWCantr...(a)sigmaxi.net> wrote:
> >
> > > On Jan 26, 4:11 am, David W. Cantrell <DWCantr...(a)sigmaxi.net> wrote:
> > > > David W.Cantrell<DWCantr...(a)sigmaxi.net> wrote:
> >
> > > > > Also see the OEIS entry for the number of decimal digits in n!:
> > > > > <http://www.research.att.com/~njas/sequences/A034886>,
> > > > > in particular, the last item in the formula section. I wonder if
> > > > > the formula given in that item is correct.
> >
> > > The author of [1] recently informed me that, as indicated in the OEIS
> > > entry, it is merely based on Stirling's approximation, without proof
> > > that it is correct for all n >= 2. But see below...
> >
> > > > Now that I've had some time to think about it, I must say that I
> > > > doubt the correctness of that formula:
> >
> > > > floor((log(2 pi n)/2 + n (log(n) - 1))/log(10)) + 1 [1]
> >
> > > > Indeed, I now have a heuristic argument that [1] should fail to
> > > > give the number of decimal digits in n! correctly for infinitely
> > > > many n, with the first failure occurring before n = 600000 with
> > > > 50% probability and before n = 6 * 10^11 with 100% probability.
> > > > (Of course, the first n for which [1] fails may be difficult to
> > > > find.)
> >
> > > My heuristic argument is very simple. Note that the number of decimal
> > > digits in n! is (obviously) given by
> > > floor(logGamma(n + 1)/log(10)) + 1. The difference between
> > > logGamma(n + 1)/log(10) and (log(2 pi n)/2 + n (log(n) - 1))/log(10)
> > > is approximately 1/(12 log(10) n), the infinite sum of which goes to
> > > infinity, albeit very slowly. (Harmonic series.) That is, the sum of
> > > the lengths of the intervals
> > > [(log(2 pi n)/2 + n (log(n) - 1))/log(10), logGamma(n + 1)/log(10)]
> > > is infinite. If any such interval happens to contain an integer, then
> > > floor of (log(2 pi n)/2 + n (log(n) - 1))/log(10) will be smaller by
> > > 1 than floor of logGamma(n + 1)/log(10), and so cause formula [1] to
> > > fail.
> >
> > > But I have not found any n for which [1] fails. And so I am now led
> > > to wonder if there is something happening which is not captured in
> > > the heuristic argument -- number theoretically, maybe -- which causes
> > > [1] to give the correct number of decial digits in n! for all n >= 2
> > > I would appreciate other people's thoughts on this matter.
> >
> > > > But there is a very simple modification to the formula which
> > > > should, according to a heuristic argument, make it valid, with high
> > > > probability, for all n >= 1:
> >
> > > > floor((log(2 pi n)/2 + n (log(n) - 1) + 1/(12 n))/log(10)) + 1 [2]
> >
> > > Using the same sort of heuristic reasoning:
> > > The difference between
> > > (log(2 pi n)/2 + n (log(n) - 1) + 1/(12 n))/log(10) and
> > > logGamma(n + 1)/log(10) is approximately 1/(360 log(10) n^3).
> > > The sum of the lengths of the intervals
> > > [logGamma(n + 1)/log(10), (log(2 pi n)/2 + n (log(n) - 1) + 1/(12
> > > n))/log(10)], from n = 1..oo, is zeta(3)/(360 log(10)) = 0.00145...,
> > > making it unlikely that any such interval would happen to contain an
> > > integer.
> >
> > > Investigations in other bases:
> >
> > > Of course, my formula [2] can be modified to give the number of
> > > base-b digits in n! by simply replacing log(10) with log(b). It seems
> > > that, so modified, [2] works correctly for all n >= 1 and for all
> > > bases b.
> >
> > > But, when similarly modified, formula [1] fails in some bases. (It
> > > fails in all bases for n = 1, but that is of little interest.
> > > Henceforth, we consider just n > 1.) It will fail at n in base b if
> > > b = n! For example, it fails in base 6 for n = 3 because b = 6 = 3!
> > > and it fails in base 2 for n = 2 because b = 2 = 2! Here are the
> > > other failures I've found in various bases:
> >
> > > b n
> > > 5 47
> > > 6 2751
> > > 12 516
> > > 14 29
> > > 14 4799
> > > 14 11321
> >
> > Oops. To avoid confusion I should have specified that my values of n
> > are given in _decimal_ notation.
> >
> > > In all those cases, not surprisingly, n! in base b consists of 1
> > > followed by some 0's. For example, in base 14, we have 29! =
> > > 1006bb5d6db12a825d9398ac0000
> >
> > That is, (29_10)! = 1006bb5d6db12a825d9398ac0000_14.
> >
> > DWC
> >
> > > The fact that formula [1], when modified for base b, fails in some
> > > bases strengthens my skepticism about its working in base 10 for
> > > all n >= 2.
>
> Hi, David:
>
> I suspect your heuristic arguments are correct,
> and we are simply frustrated by the horribly
> slow divergence of the harmonic series.

Thanks, Chip. That was my initial thinking too.

> Have the intervals of disagreement between [1]
> and [2] been exhaustively checked up to some
> point, e.g. n = 600000 as you supposed would
> give a 50% chance of counterexample? It seems
> this would be easy to program, and if none is
> found, the data might suggest some systematic
> reason for failure.

Yes, it's easy to program. Using Mathematica, apparently [1] works
correctly up to at least 10^7.

[And BTW, for the nondecimal bases, my investigations were exhaustive
up to n = 10^6 for binary through hexadecimal.]

Regards,
David
From: David W. Cantrell on 6 Feb 2010 12:45
richard(a)cogsci.ed.ac.uk (Richard Tobin) wrote:
> In article <20100206120206.850$lA(a)newsreader.com>,
> David W. Cantrell <DWCantrell(a)sigmaxi.net> wrote:
>
> >> Could you determine whether it fails for 3121515!. I don't have the
> >> facilities to hand to check it accurately.
>
> That was just the first case where 64-bit IEEE floating point numbers
> give different results, calculating the factorial by adding the logs.
> Of course, I couldn't tell whether it was the accumulated floating
> point error or a real failure.
>
> >Any other suggestions?
>
> Using "long double", the first candidates are 48655817 and 91948153
> (beyond that they occur often, suggesting that the floating point
> error has become significant).

Thanks again. The formula in OEIS also works correctly for both of those
candidates.

David
From: David W. Cantrell on 6 Feb 2010 16:52
David W. Cantrell <DWCantrell(a)sigmaxi.net> wrote:
> Chip Eastham <hardmath(a)gmail.com> wrote:
> > On Feb 6, 1:56=A0am, David W. Cantrell <DWCantr...(a)sigmaxi.net> wrote:
> > > "David W. Cantrell" <DWCantr...(a)sigmaxi.net> wrote:
> > >
> > > > On Jan 26, 4:11 am, David W. Cantrell <DWCantr...(a)sigmaxi.net>
> > > > wrote:
> > > > > David W.Cantrell<DWCantr...(a)sigmaxi.net> wrote:
> > >
> > > > > > Also see the OEIS entry for the number of decimal digits in n!:
> > > > > > <http://www.research.att.com/~njas/sequences/A034886>,
> > > > > > in particular, the last item in the formula section. I wonder
> > > > > > if the formula given in that item is correct.
> > >
> > > > The author of [1] recently informed me that, as indicated in the
> > > > OEIS entry, it is merely based on Stirling's approximation, without
> > > > proof that it is correct for all n >= 2. But see below...
> > >
> > > > > Now that I've had some time to think about it, I must say that I
> > > > > doubt the correctness of that formula:
> > >
> > > > > floor((log(2 pi n)/2 + n (log(n) - 1))/log(10)) + 1 [1]
> > >
> > > > > Indeed, I now have a heuristic argument that [1] should fail to
> > > > > give the number of decimal digits in n! correctly for infinitely
> > > > > many n, with the first failure occurring before n = 600000 with
> > > > > 50% probability and before n = 6 * 10^11 with 100% probability.
> > > > > (Of course, the first n for which [1] fails may be difficult to
> > > > > find.)
> > >
> > > > My heuristic argument is very simple. Note that the number of
> > > > decimal digits in n! is (obviously) given by
> > > > floor(logGamma(n + 1)/log(10)) + 1. The difference between
> > > > logGamma(n + 1)/log(10) and (log(2 pi n)/2 + n (log(n) -
> > > > 1))/log(10) is approximately 1/(12 log(10) n), the infinite sum of
> > > > which goes to infinity, albeit very slowly. (Harmonic series.) That
> > > > is, the sum of the lengths of the intervals
> > > > [(log(2 pi n)/2 + n (log(n) - 1))/log(10), logGamma(n + 1)/log(10)]
> > > > is infinite. If any such interval happens to contain an integer,
> > > > then floor of (log(2 pi n)/2 + n (log(n) - 1))/log(10) will be
> > > > smaller by 1 than floor of logGamma(n + 1)/log(10), and so cause
> > > > formula [1] to fail.
> > >
> > > > But I have not found any n for which [1] fails. And so I am now led
> > > > to wonder if there is something happening which is not captured in
> > > > the heuristic argument -- number theoretically, maybe -- which
> > > > causes [1] to give the correct number of decial digits in n! for
> > > > all n >= 2 I would appreciate other people's thoughts on this
> > > > matter.
> > >
> > > > > But there is a very simple modification to the formula which
> > > > > should, according to a heuristic argument, make it valid, with
> > > > > high probability, for all n >= 1:
> > >
> > > > > floor((log(2 pi n)/2 + n (log(n) - 1) + 1/(12 n))/log(10)) + 1 [2]
> > >
> > > > Using the same sort of heuristic reasoning:
> > > > The difference between
> > > > (log(2 pi n)/2 + n (log(n) - 1) + 1/(12 n))/log(10) and
> > > > logGamma(n + 1)/log(10) is approximately 1/(360 log(10) n^3).
> > > > The sum of the lengths of the intervals
> > > > [logGamma(n + 1)/log(10), (log(2 pi n)/2 + n (log(n) - 1) + 1/(12
> > > > n))/log(10)], from n = 1..oo, is zeta(3)/(360 log(10)) =
> > > > 0.00145..., making it unlikely that any such interval would happen
> > > > to contain an integer.
> > >
> > > > Investigations in other bases:
> > >
> > > > Of course, my formula [2] can be modified to give the number of
> > > > base-b digits in n! by simply replacing log(10) with log(b). It
> > > > seems that, so modified, [2] works correctly for all n >= 1 and for
> > > > all bases b.
> > >
> > > > But, when similarly modified, formula [1] fails in some bases. (It
> > > > fails in all bases for n = 1, but that is of little interest.
> > > > Henceforth, we consider just n > 1.) It will fail at n in base b if
> > > > b = n! For example, it fails in base 6 for n = 3 because b = 6 = 3!
> > > > and it fails in base 2 for n = 2 because b = 2 = 2! Here are the
> > > > other failures I've found in various bases:
> > >
> > > > b n
> > > > 5 47
> > > > 6 2751
> > > > 12 516
> > > > 14 29
> > > > 14 4799
> > > > 14 11321
> > >
> > > Oops. To avoid confusion I should have specified that my values of n
> > > are given in _decimal_ notation.
> > >
> > > > In all those cases, not surprisingly, n! in base b consists of 1
> > > > followed by some 0's. For example, in base 14, we have 29! =
> > > > 1006bb5d6db12a825d9398ac0000
> > >
> > > That is, (29_10)! = 1006bb5d6db12a825d9398ac0000_14.
> > >
> > > DWC
> > >
> > > > The fact that formula [1], when modified for base b, fails in some
> > > > bases strengthens my skepticism about its working in base 10 for
> > > > all n >= 2.
> >
> > Hi, David:
> >
> > I suspect your heuristic arguments are correct,
> > and we are simply frustrated by the horribly
> > slow divergence of the harmonic series.
>
> Thanks, Chip. That was my initial thinking too.
>
> > Have the intervals of disagreement between [1]
> > and [2] been exhaustively checked up to some
> > point, e.g. n = 600000 as you supposed would
> > give a 50% chance of counterexample? It seems
> > this would be easy to program, and if none is
> > found, the data might suggest some systematic
> > reason for failure.
>
> Yes, it's easy to program. Using Mathematica, apparently [1] works
> correctly up to at least 10^7.
>
> [And BTW, for the nondecimal bases, my investigations were exhaustive
> up to n = 10^6 for binary through hexadecimal.]

Oops, again. There was an important omission from my table of failures.

In some nondecimal bases, I searched beyond n = 10^6. In particular, in
binary, I had searched up to n = 2 * 10^8, but I had somehow overlooked a
failure which Mathematica had already found at n = 25463836:

Formula [1] {with log(10) replaced by log(2), of course} predicts that
25463836! should have 589723393 bits. But it actually has 589723394 bits.
The interval
[(log(2 pi n)/2 + n (log(n) - 1))/log(2), logGamma(n + 1)/log(2)] = [589723392.9999999961..., 589723393.00000000083...] contains an integer.
Expressing 25463836! in binary, between the leading 1 and the
next 1, there are about 30 0's, if my calculation is correct.

So, in binary, [1] fails for n = 1, 2 and 25463836. And I suppose
that there are other failures beyond n = 2 * 10^8.

David
From: Chip Eastham on 8 Feb 2010 09:15
On Feb 6, 4:52 pm, David W. Cantrell <DWCantr...(a)sigmaxi.net> wrote:
> David W. Cantrell <DWCantr...(a)sigmaxi.net> wrote:
>
>
>
>
>
> > Chip Eastham <hardm...(a)gmail.com> wrote:
> > > On Feb 6, 1:56=A0am, David W. Cantrell <DWCantr...(a)sigmaxi.net> wrote:
> > > > "David W. Cantrell" <DWCantr...(a)sigmaxi.net> wrote:
>
> > > > > On Jan 26, 4:11 am, David W. Cantrell <DWCantr...(a)sigmaxi.net>
> > > > > wrote:
> > > > > > David W.Cantrell<DWCantr...(a)sigmaxi.net> wrote:
>
> > > > > > > Also see the OEIS entry for the number of decimal digits in n!:
> > > > > > > <http://www.research.att.com/~njas/sequences/A034886>,
> > > > > > > in particular, the last item in the formula section. I wonder
> > > > > > > if the formula given in that item is correct.
>
> > > > > The author of [1] recently informed me that, as indicated in the
> > > > > OEIS entry, it is merely based on Stirling's approximation, without
> > > > > proof that it is correct for all n >= 2. But see below...
>
> > > > > > Now that I've had some time to think about it, I must say that I
> > > > > > doubt the correctness of that formula:
>
> > > > > > floor((log(2 pi n)/2 + n (log(n) - 1))/log(10)) + 1          [1]
>
> > > > > > Indeed, I now have a heuristic argument that [1] should fail to
> > > > > > give the number of decimal digits in n! correctly for infinitely
> > > > > > many n, with the first failure occurring before n = 600000 with
> > > > > > 50% probability and before n = 6 * 10^11 with 100% probability.
> > > > > > (Of course, the first n for which [1] fails may be difficult to
> > > > > > find.)
>
> > > > > My heuristic argument is very simple. Note that the number of
> > > > > decimal digits in n! is (obviously) given by
> > > > > floor(logGamma(n + 1)/log(10)) + 1. The difference between
> > > > > logGamma(n + 1)/log(10) and (log(2 pi n)/2 + n (log(n) -
> > > > > 1))/log(10) is approximately 1/(12 log(10) n), the infinite sum of
> > > > > which goes to infinity, albeit very slowly. (Harmonic series.) That
> > > > > is, the sum of the lengths of the intervals
> > > > > [(log(2 pi n)/2 + n (log(n) - 1))/log(10), logGamma(n + 1)/log(10)]
> > > > > is infinite. If any such interval happens to contain an integer,
> > > > > then floor of (log(2 pi n)/2 + n (log(n) - 1))/log(10) will be
> > > > > smaller by 1 than floor of logGamma(n + 1)/log(10), and so cause
> > > > > formula [1] to fail.
>
> > > > > But I have not found any n for which [1] fails. And so I am now led
> > > > > to wonder if there is something happening which is not captured in
> > > > > the heuristic argument -- number theoretically, maybe -- which
> > > > > causes [1] to give the correct number of decial digits in n! for
> > > > > all n >= 2 I would appreciate other people's thoughts on this
> > > > > matter.
>
> > > > > > But there is a very simple modification to the formula which
> > > > > > should, according to a heuristic argument, make it valid, with
> > > > > > high probability, for all n >= 1:
>
> > > > > > floor((log(2 pi n)/2 + n (log(n) - 1) + 1/(12 n))/log(10)) + 1   [2]
>
> > > > > Using the same sort of heuristic reasoning:
> > > > > The difference between
> > > > > (log(2 pi n)/2 + n (log(n) - 1) + 1/(12 n))/log(10) and
> > > > > logGamma(n + 1)/log(10) is approximately 1/(360 log(10) n^3).
> > > > > The sum of the lengths of the intervals
> > > > > [logGamma(n + 1)/log(10), (log(2 pi n)/2 + n (log(n) - 1) + 1/(12
> > > > > n))/log(10)], from n = 1..oo, is zeta(3)/(360 log(10)) =
> > > > > 0.00145..., making it unlikely that any such interval would happen
> > > > > to contain an integer.
>
> > > > > Investigations in other bases:
>
> > > > > Of course, my formula [2] can be modified to give the number of
> > > > > base-b digits in n! by simply replacing log(10) with log(b). It
> > > > > seems that, so modified, [2] works correctly for all n >= 1 and for
> > > > > all bases b.
>
> > > > > But, when similarly modified, formula [1] fails in some bases. (It
> > > > > fails in all bases for n = 1, but that is of little interest.
> > > > > Henceforth, we consider just n > 1.) It will fail at n in base b if
> > > > > b = n! For example, it fails in base 6 for n = 3 because b = 6 = 3!
> > > > > and it fails in base 2 for n = 2 because b = 2 = 2! Here are the
> > > > > other failures I've found in various bases:
>
> > > > > b   n
> > > > > 5   47
> > > > > 6   2751
> > > > > 12  516
> > > > > 14  29
> > > > > 14  4799
> > > > > 14  11321
>
> > > > Oops. To avoid confusion I should have specified that my values of n
> > > > are given in _decimal_ notation.
>
> > > > > In all those cases, not surprisingly, n! in base b consists of 1
> > > > > followed by some 0's. For example, in base 14, we have 29! =
> > > > > 1006bb5d6db12a825d9398ac0000
>
> > > > That is, (29_10)! = 1006bb5d6db12a825d9398ac0000_14.
>
> > > > DWC
>
> > > > > The fact that formula [1], when modified for base b, fails in some
> > > > > bases strengthens my skepticism about its working in base 10 for
> > > > > all n >= 2.
>
> > > Hi, David:
>
> > > I suspect your heuristic arguments are correct,
> > > and we are simply frustrated by the horribly
> > > slow divergence of the harmonic series.
>
> > Thanks, Chip. That was my initial thinking too.
>
> > > Have the intervals of disagreement between [1]
> > > and [2] been exhaustively checked up to some
> > > point, e.g. n = 600000 as you supposed would
> > > give a 50% chance of counterexample?  It seems
> > > this would be easy to program, and if none is
> > > found, the data might suggest some systematic
> > > reason for failure.
>
> > Yes, it's easy to program. Using Mathematica, apparently [1] works
> > correctly up to at least 10^7.
>
> > [And BTW, for the nondecimal bases, my investigations were exhaustive
> > up to n = 10^6 for binary through hexadecimal.]
>
> Oops, again. There was an important omission from my table of failures.
>
> In some nondecimal bases, I searched beyond n = 10^6. In particular, in
> binary, I had searched up to n = 2 * 10^8, but I had somehow overlooked a
> failure which Mathematica had already found at n = 25463836:
>
> Formula [1]  {with log(10) replaced by log(2), of course} predicts that
> 25463836! should have 589723393 bits. But it actually has 589723394 bits.
> The interval
> [(log(2 pi n)/2 + n (log(n) - 1))/log(2), logGamma(n + 1)/log(2)] = [589723392.9999999961..., 589723393.00000000083...] contains an integer.
> Expressing 25463836! in binary, between the leading 1 and the
> next 1, there are about 30 0's, if my calculation is correct.
>
> So, in binary, [1] fails for n = 1, 2 and 25463836. And I suppose
> that there are other failures beyond n = 2 * 10^8.
>
> David

To summarize briefly, we seek integer n > 1 such that interval:

[ (log(2 pi n)/2 + n (log(n) - 1))/log(10),
(log(2 pi n)/2 + n (log(n) - 1) + 1/(12n))/log(10) ]

contains an integer. My question concerns the search method.
Are you iterating over n?

I'm curious whether defining:

f(n) = (log(2 pi n)/2 + n (log(n) - 1))/log(10)

and "solving" f(n) = M for increasing values of M will give
a more efficient search method. Am I off-base, or a day late
(and dollar short) with this idea?

regards, chip
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4
Prev: Cubic equations and companion matrices
Next: Why does JH Conway say -1 is prime?