1. Forum
  2. Usenet
  3. COMP.SOFT-SYS.MATH.MAPLE

  • Do you trust Maple's results?

    From Peter Luschny@21:1/5 to All on Sat Sep 23 06:46:34 2017
    a := n -> denom(GAMMA(n+1/2)^2/(2*n*Pi)):
    b := n -> 2^(2*n+1+padic:-ordp(n,2)):
    seq(round(evalf(a(n)/b(n),800)), n=1..290);

    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From acer@21:1/5 to Peter Luschny on Mon Sep 25 13:39:47 2017
    On Saturday, September 23, 2017 at 9:46:35 AM UTC-4, Peter Luschny wrote:
    a := n -> denom(GAMMA(n+1/2)^2/(2*n*Pi)):
    b := n -> 2^(2*n+1+padic:-ordp(n,2)):
    seq(round(evalf(a(n)/b(n),800)), n=1..290);

    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0



    For n::posint larger than 100 the call GAMMA(n) returns unevaluated. (You'd likely have discover this if you didn't wrap with `denom`, to look at the intermediate result, BTW.) The help page ?GAMMA indicates that the `expand` command can deal with that
    case. But for n=201/2 a rational that doesn't do the job.

    From Maple 18 onwards, this formula could be obtained:

    lprint( convert(GAMMA(n), doublefactorial) );
    2*doublefactorial(2*n-2)*(2/Pi)^(-1/4+1/4*cos(2*Pi*n))/(2^n)

    And that `doublefactorial` could also be "expanded", say with further conversion to `factorial`.

    ee:=GAMMA(101+1/2);
    ee := GAMMA(203/2)

    ff:=convert(convert(ee,doublefactorial),factorial):

    lprint(%);
    13399280426846637026266542163537418668282223673245521\ 36009000175903930373973956200201339808971908813302485\ 337859343555595880232757068702494068539027184400769291\ 974725364464819431304931640625/25353012004564588029934\
    06410752/(1/Pi)^(1/2)

    It's interesting that this call to `simplify` works, even if conversion of `ee` is a more obscure task.

    simplify(ee - ff);
    0

    If your Maple version is even older, and does not provide that conversion formula of GAMMA(n) to that formula involving doublefactorial, here it is:

    lprint( convert(GAMMA(n), doublefactorial) ); 2*doublefactorial(2*n-2)*(2/Pi)^(-1/4+1/4*cos(2*Pi*n))/(2^n)


    In modern Maple then, one could do this:

    a := n -> GAMMA(n+1/2)^2/(2*n*Pi):
    b := n -> 2^(2*n+1+padic:-ordp(n,2)):
    c := n -> convert(convert(GAMMA(n+1/2),doublefactorial),factorial)^2/(2*n*Pi):

    seq( simplify(a(n)-c(n)), n=98..101 );
    0, 0, 0, 0

    seq(denom(c(n))/b(n), n=98..101);
    1, 1, 1, 1

    seq(denom(c(n))/b(n), n=1..290, 10);
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Peter Luschny@21:1/5 to All on Tue Sep 26 13:59:16 2017
    For n::posint larger than 100 the call GAMMA(n) returns unevaluated.
    [...] With that in mind your use of `denom` above in procedure `a`
    presumes a particular form that you might have realized would not attain.

    Acer, I'm sure I've fallen into the same trap with Maple
    several times. Besides my bad memory, there's another reason:

    It is beyond my imagination that in 2017 a CAS that wants to
    be more than a student project has such ridiculous limitations.

    And yes, this may be pointed out in the Maple documentation,
    but I don't feel like looking into things that can be taken
    for granted, if they might have an exception in Maple.

    Maybe I'm just spoiled by SageMath:

    def a(n): return (gamma(n+1/2)^2/(2*n*pi)).denominator()
    def b(n): return 2^(2*n+1+n.valuation(2))
    print [a(n)/b(n) for n in (1..290)]

    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From acer@21:1/5 to acer on Tue Sep 26 11:12:26 2017
    On Monday, September 25, 2017 at 4:39:48 PM UTC-4, acer wrote:
    On Saturday, September 23, 2017 at 9:46:35 AM UTC-4, Peter Luschny wrote:
    a := n -> denom(GAMMA(n+1/2)^2/(2*n*Pi)):
    b := n -> 2^(2*n+1+padic:-ordp(n,2)):
    seq(round(evalf(a(n)/b(n),800)), n=1..290);

    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0



    For n::posint larger than 100 the call GAMMA(n) returns unevaluated. (You'd likely have discover this if you didn't wrap with `denom`, to look at the intermediate result, BTW.) The help page ?GAMMA indicates that the `expand` command can deal with that
    case. But for n=201/2 a rational that doesn't do the job.

    From Maple 18 onwards, this formula could be obtained:

    lprint( convert(GAMMA(n), doublefactorial) );
    2*doublefactorial(2*n-2)*(2/Pi)^(-1/4+1/4*cos(2*Pi*n))/(2^n)

    And that `doublefactorial` could also be "expanded", say with further conversion to `factorial`.

    ee:=GAMMA(101+1/2);
    ee := GAMMA(203/2)

    ff:=convert(convert(ee,doublefactorial),factorial):

    lprint(%);
    13399280426846637026266542163537418668282223673245521\ 36009000175903930373973956200201339808971908813302485\ 337859343555595880232757068702494068539027184400769291\ 974725364464819431304931640625/25353012004564588029934\
    06410752/(1/Pi)^(1/2)

    It's interesting that this call to `simplify` works, even if conversion of `ee` is a more obscure task.

    simplify(ee - ff);
    0

    If your Maple version is even older, and does not provide that conversion formula of GAMMA(n) to that formula involving doublefactorial, here it is:

    lprint( convert(GAMMA(n), doublefactorial) ); 2*doublefactorial(2*n-2)*(2/Pi)^(-1/4+1/4*cos(2*Pi*n))/(2^n)


    In modern Maple then, one could do this:

    a := n -> GAMMA(n+1/2)^2/(2*n*Pi):
    b := n -> 2^(2*n+1+padic:-ordp(n,2)):
    c := n -> convert(convert(GAMMA(n+1/2),doublefactorial),factorial)^2/(2*n*Pi):

    seq( simplify(a(n)-c(n)), n=98..101 );
    0, 0, 0, 0

    seq(denom(c(n))/b(n), n=98..101);
    1, 1, 1, 1

    seq(denom(c(n))/b(n), n=1..290, 10);
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1


    I looked at some older posts in this forum and found this:
    https://groups.google.com/forum/#!topic/comp.soft-sys.math.maple/DbHyoWjEXfw

    So it seems that at least by Nov 2015 you were aware that GAMMA(201/2) would return unevaluated. With that in mind your use of `denom` above in procedure `a` presumes a particular form that you might have realized would not attain.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)