Computing the twin prime constant with Mathematica

The name of the pictureThe name of the pictureThe name of the pictureClash Royale CLAN TAG#URR8PPP











up vote
1
down vote

favorite
1












The twin prime constant is defined as



$$
C _2 = prod limits_p geqslant 3 , fracp(p-2)(p-1)^2.
$$



Translating this directly in Mathematica code gives:



N[Product[Prime[i] (Prime[i] - 2)/(Prime[i] - 1)^2, i, 2, Infinity]]


However, when I execute this code I get the following error messages:



enter image description here



What am I doing wrong?






code-review prime-numbers







share|improve this question

















  • 2




    It's just that some of the algorithms used internally by NProduct sometimes try to evaluate terms for inexact arguments, so a more clever approach is necessary.
    – J. M. is computer-less♦
    2 hours ago






  • 1




    I didn't get this message. However I suggest not using Infinity inside Product, as there will possibly be no results. Try a big number instead, like 100000
    – t-smart
    1 hour ago






  • 1




    This evaluates fine if you simply use a really big number, like 10^6
    – Carl Lange
    1 hour ago






  • 3




    See oeis.org/A005597 for better methods, how to compute this constant.
    – Vaclav Kotesovec
    1 hour ago














up vote
1
down vote

favorite
1












The twin prime constant is defined as



$$
C _2 = prod limits_p geqslant 3 , fracp(p-2)(p-1)^2.
$$



Translating this directly in Mathematica code gives:



N[Product[Prime[i] (Prime[i] - 2)/(Prime[i] - 1)^2, i, 2, Infinity]]


However, when I execute this code I get the following error messages:



enter image description here



What am I doing wrong?






code-review prime-numbers







share|improve this question

















  • 2




    It's just that some of the algorithms used internally by NProduct sometimes try to evaluate terms for inexact arguments, so a more clever approach is necessary.
    – J. M. is computer-less♦
    2 hours ago






  • 1




    I didn't get this message. However I suggest not using Infinity inside Product, as there will possibly be no results. Try a big number instead, like 100000
    – t-smart
    1 hour ago






  • 1




    This evaluates fine if you simply use a really big number, like 10^6
    – Carl Lange
    1 hour ago






  • 3




    See oeis.org/A005597 for better methods, how to compute this constant.
    – Vaclav Kotesovec
    1 hour ago












up vote
1
down vote

favorite
1









up vote
1
down vote

favorite
1






1





The twin prime constant is defined as



$$
C _2 = prod limits_p geqslant 3 , fracp(p-2)(p-1)^2.
$$



Translating this directly in Mathematica code gives:



N[Product[Prime[i] (Prime[i] - 2)/(Prime[i] - 1)^2, i, 2, Infinity]]


However, when I execute this code I get the following error messages:



enter image description here



What am I doing wrong?






code-review prime-numbers







share|improve this question













The twin prime constant is defined as



$$
C _2 = prod limits_p geqslant 3 , fracp(p-2)(p-1)^2.
$$



Translating this directly in Mathematica code gives:



N[Product[Prime[i] (Prime[i] - 2)/(Prime[i] - 1)^2, i, 2, Infinity]]


However, when I execute this code I get the following error messages:



enter image description here



What am I doing wrong?






code-review prime-numbers




code-review prime-numbers






share|improve this question













share|improve this question











share|improve this question




share|improve this question










asked 2 hours ago









PierreTheFermented

21717




21717







  • 2




    It's just that some of the algorithms used internally by NProduct sometimes try to evaluate terms for inexact arguments, so a more clever approach is necessary.
    – J. M. is computer-less♦
    2 hours ago






  • 1




    I didn't get this message. However I suggest not using Infinity inside Product, as there will possibly be no results. Try a big number instead, like 100000
    – t-smart
    1 hour ago






  • 1




    This evaluates fine if you simply use a really big number, like 10^6
    – Carl Lange
    1 hour ago






  • 3




    See oeis.org/A005597 for better methods, how to compute this constant.
    – Vaclav Kotesovec
    1 hour ago












  • 2




    It's just that some of the algorithms used internally by NProduct sometimes try to evaluate terms for inexact arguments, so a more clever approach is necessary.
    – J. M. is computer-less♦
    2 hours ago






  • 1




    I didn't get this message. However I suggest not using Infinity inside Product, as there will possibly be no results. Try a big number instead, like 100000
    – t-smart
    1 hour ago






  • 1




    This evaluates fine if you simply use a really big number, like 10^6
    – Carl Lange
    1 hour ago






  • 3




    See oeis.org/A005597 for better methods, how to compute this constant.
    – Vaclav Kotesovec
    1 hour ago







2




2




It's just that some of the algorithms used internally by NProduct sometimes try to evaluate terms for inexact arguments, so a more clever approach is necessary.
– J. M. is computer-less♦
2 hours ago




It's just that some of the algorithms used internally by NProduct sometimes try to evaluate terms for inexact arguments, so a more clever approach is necessary.
– J. M. is computer-less♦
2 hours ago




1




1




I didn't get this message. However I suggest not using Infinity inside Product, as there will possibly be no results. Try a big number instead, like 100000
– t-smart
1 hour ago




I didn't get this message. However I suggest not using Infinity inside Product, as there will possibly be no results. Try a big number instead, like 100000
– t-smart
1 hour ago




1




1




This evaluates fine if you simply use a really big number, like 10^6
– Carl Lange
1 hour ago




This evaluates fine if you simply use a really big number, like 10^6
– Carl Lange
1 hour ago




3




3




See oeis.org/A005597 for better methods, how to compute this constant.
– Vaclav Kotesovec
1 hour ago




See oeis.org/A005597 for better methods, how to compute this constant.
– Vaclav Kotesovec
1 hour ago










2 Answers
2






active

oldest

votes

















up vote
3
down vote



accepted










The problem with your code is that N[Product[...]] calls NProduct[...] which turns integers into reals before Prime can evaluate, causing the error: See Details and Options section of documentation of NProduct.



For your calculation, you actually do not need these. It looks like the result converges around $0.660162$:



Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k]], k, 5]



0.665138, 0.66033, 0.66017, 0.660162, 0.660162




Since the product is monotonic and you are using N at the end, it is not necessary to carry out whole multiplication all the way upto infinity: Depending on the interested accuracy, you can go to higher orders as well. For example,



Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k], 8], k, 5]


yields




0.66513840, 0.66033029, 0.66017020, 0.66016232, 0.66016185




which reveals that using only first $10^5$ primes, which is all primes upto $1299709$, is sufficient to get a result with error of the order $10^-6$.



For the record: Mathematica fails to do the calculation analytically:



Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, Infinity]



$prod _i=2^infty fracleft(p_i-2right) p_ileft(p_i-1right)^2$







share|improve this answer




















  • Thank you. So if I want an error of maximum $10^-x$, I only need to use the first $10^x-1$ primes?
    – PierreTheFermented
    25 mins ago










  • @PierreTheFermented No, not necessarily. I was just trying to say that you do not need all primes and depending on the accuracy needed, different number of primes will be sufficient. I am not sure if there is an analytical way to predict that number though. In any case, using more than first 10^6 primes will probably be fairly slow.
    – Soner
    20 mins ago


















up vote
3
down vote













As not all primes are known, Mathematica will have a problem to compute the constant exactly. But the sequence appears to converge quickly to some number $0.660162...$. (Soner posted that already.)



So let me give a way to compute the sequence of partial products in a more efficient way.



n = 10000;
seq = FoldList[#1 #2 (#2 - 2)/(#2 - 1)^2 &, 1, Prime[Range[2, n]]];
N[seq[[-1]]]
ListLinePlot[N[seq], ConstantArray[N[seq[[-1]]], Length[seq]],
 PlotRange -> 1/2, 1,
 PlotStyle -> Thick, Dashed
 ]


enter image description here






share|improve this answer




















    Your Answer




    StackExchange.ifUsing("editor", function ()
    return StackExchange.using("mathjaxEditing", function ()
    StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix)
    StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
    );
    );
    , "mathjax-editing");

    StackExchange.ready(function()
    var channelOptions =
    tags: "".split(" "),
    id: "387"
    ;
    initTagRenderer("".split(" "), "".split(" "), channelOptions);

    StackExchange.using("externalEditor", function()
    // Have to fire editor after snippets, if snippets enabled
    if (StackExchange.settings.snippets.snippetsEnabled)
    StackExchange.using("snippets", function()
    createEditor();
    );

    else
    createEditor();

    );

    function createEditor()
    StackExchange.prepareEditor(
    heartbeatType: 'answer',
    convertImagesToLinks: false,
    noModals: false,
    showLowRepImageUploadWarning: true,
    reputationToPostImages: null,
    bindNavPrevention: true,
    postfix: "",
    onDemand: true,
    discardSelector: ".discard-answer"
    ,immediatelyShowMarkdownHelp:true
    );



    );













     

    draft saved


    draft discarded


















    StackExchange.ready(
    function ()
    StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f184324%2fcomputing-the-twin-prime-constant-with-mathematica%23new-answer', 'question_page');

    );

    Post as a guest

















    2 Answers
    2






    active

    oldest

    votes








    2 Answers
    2






    active

    oldest

    votes









    active

    oldest

    votes






    active

    oldest

    votes








    up vote
    3
    down vote



    accepted










    The problem with your code is that N[Product[...]] calls NProduct[...] which turns integers into reals before Prime can evaluate, causing the error: See Details and Options section of documentation of NProduct.



    For your calculation, you actually do not need these. It looks like the result converges around $0.660162$:



    Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k]], k, 5]



    0.665138, 0.66033, 0.66017, 0.660162, 0.660162




    Since the product is monotonic and you are using N at the end, it is not necessary to carry out whole multiplication all the way upto infinity: Depending on the interested accuracy, you can go to higher orders as well. For example,



    Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k], 8], k, 5]


    yields




    0.66513840, 0.66033029, 0.66017020, 0.66016232, 0.66016185




    which reveals that using only first $10^5$ primes, which is all primes upto $1299709$, is sufficient to get a result with error of the order $10^-6$.



    For the record: Mathematica fails to do the calculation analytically:



    Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, Infinity]



    $prod _i=2^infty fracleft(p_i-2right) p_ileft(p_i-1right)^2$







    share|improve this answer




















    • Thank you. So if I want an error of maximum $10^-x$, I only need to use the first $10^x-1$ primes?
      – PierreTheFermented
      25 mins ago










    • @PierreTheFermented No, not necessarily. I was just trying to say that you do not need all primes and depending on the accuracy needed, different number of primes will be sufficient. I am not sure if there is an analytical way to predict that number though. In any case, using more than first 10^6 primes will probably be fairly slow.
      – Soner
      20 mins ago















    up vote
    3
    down vote



    accepted










    The problem with your code is that N[Product[...]] calls NProduct[...] which turns integers into reals before Prime can evaluate, causing the error: See Details and Options section of documentation of NProduct.



    For your calculation, you actually do not need these. It looks like the result converges around $0.660162$:



    Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k]], k, 5]



    0.665138, 0.66033, 0.66017, 0.660162, 0.660162




    Since the product is monotonic and you are using N at the end, it is not necessary to carry out whole multiplication all the way upto infinity: Depending on the interested accuracy, you can go to higher orders as well. For example,



    Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k], 8], k, 5]


    yields




    0.66513840, 0.66033029, 0.66017020, 0.66016232, 0.66016185




    which reveals that using only first $10^5$ primes, which is all primes upto $1299709$, is sufficient to get a result with error of the order $10^-6$.



    For the record: Mathematica fails to do the calculation analytically:



    Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, Infinity]



    $prod _i=2^infty fracleft(p_i-2right) p_ileft(p_i-1right)^2$







    share|improve this answer




















    • Thank you. So if I want an error of maximum $10^-x$, I only need to use the first $10^x-1$ primes?
      – PierreTheFermented
      25 mins ago










    • @PierreTheFermented No, not necessarily. I was just trying to say that you do not need all primes and depending on the accuracy needed, different number of primes will be sufficient. I am not sure if there is an analytical way to predict that number though. In any case, using more than first 10^6 primes will probably be fairly slow.
      – Soner
      20 mins ago













    up vote
    3
    down vote



    accepted







    up vote
    3
    down vote



    accepted






    The problem with your code is that N[Product[...]] calls NProduct[...] which turns integers into reals before Prime can evaluate, causing the error: See Details and Options section of documentation of NProduct.



    For your calculation, you actually do not need these. It looks like the result converges around $0.660162$:



    Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k]], k, 5]



    0.665138, 0.66033, 0.66017, 0.660162, 0.660162




    Since the product is monotonic and you are using N at the end, it is not necessary to carry out whole multiplication all the way upto infinity: Depending on the interested accuracy, you can go to higher orders as well. For example,



    Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k], 8], k, 5]


    yields




    0.66513840, 0.66033029, 0.66017020, 0.66016232, 0.66016185




    which reveals that using only first $10^5$ primes, which is all primes upto $1299709$, is sufficient to get a result with error of the order $10^-6$.



    For the record: Mathematica fails to do the calculation analytically:



    Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, Infinity]



    $prod _i=2^infty fracleft(p_i-2right) p_ileft(p_i-1right)^2$







    share|improve this answer












    The problem with your code is that N[Product[...]] calls NProduct[...] which turns integers into reals before Prime can evaluate, causing the error: See Details and Options section of documentation of NProduct.



    For your calculation, you actually do not need these. It looks like the result converges around $0.660162$:



    Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k]], k, 5]



    0.665138, 0.66033, 0.66017, 0.660162, 0.660162




    Since the product is monotonic and you are using N at the end, it is not necessary to carry out whole multiplication all the way upto infinity: Depending on the interested accuracy, you can go to higher orders as well. For example,



    Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, 10^k], 8], k, 5]


    yields




    0.66513840, 0.66033029, 0.66017020, 0.66016232, 0.66016185




    which reveals that using only first $10^5$ primes, which is all primes upto $1299709$, is sufficient to get a result with error of the order $10^-6$.



    For the record: Mathematica fails to do the calculation analytically:



    Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, i, 2, Infinity]



    $prod _i=2^infty fracleft(p_i-2right) p_ileft(p_i-1right)^2$








    share|improve this answer












    share|improve this answer



    share|improve this answer










    answered 1 hour ago









    Soner

    62839




    62839











    • Thank you. So if I want an error of maximum $10^-x$, I only need to use the first $10^x-1$ primes?
      – PierreTheFermented
      25 mins ago










    • @PierreTheFermented No, not necessarily. I was just trying to say that you do not need all primes and depending on the accuracy needed, different number of primes will be sufficient. I am not sure if there is an analytical way to predict that number though. In any case, using more than first 10^6 primes will probably be fairly slow.
      – Soner
      20 mins ago

















    • Thank you. So if I want an error of maximum $10^-x$, I only need to use the first $10^x-1$ primes?
      – PierreTheFermented
      25 mins ago










    • @PierreTheFermented No, not necessarily. I was just trying to say that you do not need all primes and depending on the accuracy needed, different number of primes will be sufficient. I am not sure if there is an analytical way to predict that number though. In any case, using more than first 10^6 primes will probably be fairly slow.
      – Soner
      20 mins ago
















    Thank you. So if I want an error of maximum $10^-x$, I only need to use the first $10^x-1$ primes?
    – PierreTheFermented
    25 mins ago




    Thank you. So if I want an error of maximum $10^-x$, I only need to use the first $10^x-1$ primes?
    – PierreTheFermented
    25 mins ago












    @PierreTheFermented No, not necessarily. I was just trying to say that you do not need all primes and depending on the accuracy needed, different number of primes will be sufficient. I am not sure if there is an analytical way to predict that number though. In any case, using more than first 10^6 primes will probably be fairly slow.
    – Soner
    20 mins ago





    @PierreTheFermented No, not necessarily. I was just trying to say that you do not need all primes and depending on the accuracy needed, different number of primes will be sufficient. I am not sure if there is an analytical way to predict that number though. In any case, using more than first 10^6 primes will probably be fairly slow.
    – Soner
    20 mins ago











    up vote
    3
    down vote













    As not all primes are known, Mathematica will have a problem to compute the constant exactly. But the sequence appears to converge quickly to some number $0.660162...$. (Soner posted that already.)



    So let me give a way to compute the sequence of partial products in a more efficient way.



    n = 10000;
    seq = FoldList[#1 #2 (#2 - 2)/(#2 - 1)^2 &, 1, Prime[Range[2, n]]];
    N[seq[[-1]]]
    ListLinePlot[N[seq], ConstantArray[N[seq[[-1]]], Length[seq]],
     PlotRange -> 1/2, 1,
     PlotStyle -> Thick, Dashed
     ]


    enter image description here






    share|improve this answer
























      up vote
      3
      down vote













      As not all primes are known, Mathematica will have a problem to compute the constant exactly. But the sequence appears to converge quickly to some number $0.660162...$. (Soner posted that already.)



      So let me give a way to compute the sequence of partial products in a more efficient way.



      n = 10000;
      seq = FoldList[#1 #2 (#2 - 2)/(#2 - 1)^2 &, 1, Prime[Range[2, n]]];
      N[seq[[-1]]]
      ListLinePlot[N[seq], ConstantArray[N[seq[[-1]]], Length[seq]],
       PlotRange -> 1/2, 1,
       PlotStyle -> Thick, Dashed
       ]


      enter image description here






      share|improve this answer






















        up vote
        3
        down vote










        up vote
        3
        down vote









        As not all primes are known, Mathematica will have a problem to compute the constant exactly. But the sequence appears to converge quickly to some number $0.660162...$. (Soner posted that already.)



        So let me give a way to compute the sequence of partial products in a more efficient way.



        n = 10000;
        seq = FoldList[#1 #2 (#2 - 2)/(#2 - 1)^2 &, 1, Prime[Range[2, n]]];
        N[seq[[-1]]]
        ListLinePlot[N[seq], ConstantArray[N[seq[[-1]]], Length[seq]],
         PlotRange -> 1/2, 1,
         PlotStyle -> Thick, Dashed
         ]


        enter image description here






        share|improve this answer












        As not all primes are known, Mathematica will have a problem to compute the constant exactly. But the sequence appears to converge quickly to some number $0.660162...$. (Soner posted that already.)



        So let me give a way to compute the sequence of partial products in a more efficient way.



        n = 10000;
        seq = FoldList[#1 #2 (#2 - 2)/(#2 - 1)^2 &, 1, Prime[Range[2, n]]];
        N[seq[[-1]]]
        ListLinePlot[N[seq], ConstantArray[N[seq[[-1]]], Length[seq]],
         PlotRange -> 1/2, 1,
         PlotStyle -> Thick, Dashed
         ]


        enter image description here







        share|improve this answer












        share|improve this answer



        share|improve this answer










        answered 1 hour ago









        Henrik Schumacher

        42.5k261125




        42.5k261125



























             

            draft saved


            draft discarded















































             


            draft saved


            draft discarded














            StackExchange.ready(
            function ()
            StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f184324%2fcomputing-the-twin-prime-constant-with-mathematica%23new-answer', 'question_page');

            );

            Post as a guest
































































            Comments

            Post a Comment

            Popular posts from this blog

            What does second last employer means? [closed]

            Image
            Clash Royale CLAN TAG #URR8PPP .everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty margin-bottom:0; up vote -3 down vote favorite Say, for example I have worked in the following companies, Company A => 2010 - 2011 Company B => 2011 - 2012 Company C => 2012 - Till Date Company D => Have offer Now assume Company D is asking to submit documents from my second last employer. So from the list of companies given, which company would be qualified as second last employer? software-industry employer share | improve this question asked Aug 10 '14 at 5:16 Rajaprabhu Aravindasamy 559 1 6 13 closed as off-topic by Jim G., jmort253 ♦ Aug 11 '14 at 1:58 This question appears to be off-topic. The users who voted to close gave this specific reason: "Questions seeking advice on company-specific regulations, agreements, or policies should be directed to your manager or HR department. Questions that address...
            Read more

            List of Gilmore Girls characters

            Image
            The season three cast—(Back row): Lane, Michel, Paris, Emily, Richard, Sookie, Miss Patty, Kirk; Front row: Jess, Luke, Lorelai, Rory, Dean This is a list of characters for the comedy-drama television series Gilmore Girls . Contents 1 Main characters 1.1 Lorelai Gilmore 1.2 Rory Gilmore 1.3 Sookie St. James 1.4 Lane Kim 1.5 Michel Gerard 1.6 Luke Danes 1.7 Emily Gilmore 1.8 Richard Gilmore 1.9 Paris Geller 1.10 Dean Forester 1.11 Jess Mariano 1.12 Kirk Gleason 1.13 Jason Stiles 1.14 Logan Huntzberger 2 Major recurring characters 2.1 Christopher Hayden 2.2 Jackson Belleville 2.3 Mrs. Kim 2.4 Tristin DuGray 2.5 Max Medina 2.6 Taylor Doose 2.7 Dave Rygalski 2.8 Marty 3 Recurring characters 3.1 Stars Hollow 3.2 Chilton 3.3 Yale 3.4 Other 4 Notable guest stars 5 Unseen characters 6 Notes 7 References Main characters Actor Character Appearances Season 1 Season 2 Season 3 Season 4 Season 5 Season 6 Season 7 A Year in the Life Lauren Graham Lorelai Gilmore...
            Read more

            Confectionery

            Image
            Prepared foods made with a lot of sugar or other cabohydrates "Sweetmeat" redirects here. For the racehorse, see Sweetmeat (horse). Not to be confused with Sweetbread. This Kransekake is a traditional Scandinavian baker's confection. Confectionery is the art of making confections , which are food items that are rich in sugar and carbohydrates. Exact definitions are difficult. [1] In general, though, confectionery is divided into two broad and somewhat overlapping categories, bakers' confections and sugar confections. [2] Bakers' confectionery, also called flour confections , includes principally sweet pastries, cakes, and similar baked goods. Sugar confectionery includes candies ( sweets in British English), candied nuts, chocolates, chewing gum, bubble gum, pastillage, and other confections that are made primarily of sugar. In some cases, chocolate confections (confections made of chocolate) are treated as a separate category, as are sugar-free versions of ...
            Read more