The previous timespacex comparison may not seem fair, because the q: method takes 1e5 ints, while the zeta method takes 100, but that's because convergence is very slow for the q: method mean >./@(0&,)@(_&q:)"0 }.i.100 1.63636
First 100 ints, gives a value nowhere near nivens constnat mean >./@(0&,)@(_&q:)"0 }.i.1e6 1.70519 The first 1 million ints gives a closer approximation. -------------------------------------------- On Wed, 7/6/16, 'Jon Hough' via Programming <[email protected]> wrote: Subject: Re: [Jprogramming] Niven's constant To: [email protected] Date: Wednesday, July 6, 2016, 8:37 PM I'm sure this isn't what you are looking for but I had some time today, so I calculated Nivens constant using the zeta function and bernoullis. My guess is that it is computationally less expensive than using q: Step 1: Generate Bernoulli Numbers. We get a list of (even) Bernoulli numbers up to the given number. NB. Generates all even bernoulli numbers up to y generatebernoulli=: 3 : 0 arr=. (y + 2) $ 0 arr=. (1) 1} arr bernoulliList=. 1 flag=. 1 NB. true a=. 1 b=. 1 c=. _2 d=. 1 size=. (2 * y) - 2 for_j. i. size do. if. flag = 1 do. a=. a + 1 b=. 4 * b c=. _1 * c d=. c * (b - 1) k=. a while. k > 0 do. tmp0=. k { arr tmp1=. (k+1) { arr arr=. (tmp0 + tmp1) k} arr k=. k - 1 end. else. k=. 1 while. k < a do. tmp0=. k { arr tmp1=. (k-1) { arr arr=. (tmp0 + tmp1) k} arr k=. >: k end. bernoulliList=. bernoulliList, (x: ((a-1){arr) % d) end. flag=. 2 | >: flag end. bernoulliList ) STEP 2: Use the relation between bernoulli numbers and zeta to write a verb that computes the zeta (2 *y), for argument y. left argument, x, should be the previously calculated bernoulli number list. NB. Calculates Even Zeta, i.e. Zeta (2 *y). NB. Uses the relation between even zeta values and bernoulli numbers. calcevenzeta =: 4 : 0 bernoullis =. x m =. 2 * y (_1^(y+1)) * ((y{bernoullis) * (2 * o. 1)^m) % (2 * ! m) ) STEP 3: For odd numbers we have to do raw computation. So compute zeta(y) for 200 terms. NB. Calculates zeta (y), but should be used for odd y calcoddzeta =: 3 : 0 arr =.%>: i. 200 NB. 200 is enough +/arr ^ y ) STEP 4: calculate Nivens constant bernoullilist =: generatebernoulli 100 NB. we don't even need this many zetaapproximations =:(bl&calcevenzeta@:-:)`calcoddzeta@.(1&=(2&|)) NB. two cases, even or odd integers. NIVENS_CONSTANT =: 1+ +/ 1- % zetaapproximations"(0) 2 + i. 100 NB. we will approximate with the first 100 integers, but I doubt this many is necessary. this gives 1.7052 I wrapped all the above in a single verb and ran timespacex timespacex 'allcalc 100' NB. calculate first 100 terms only. 0.044685 34240 Doing raw calcs with q: is way more expensive: timespacex 'mean >./@(0&,)@(_&q:)"0 }.i.1e5' 1.5124 9.04083e6 -------------------------------------------- On Wed, 7/6/16, jonghough via Programming <[email protected]> wrote: Subject: Re: [Jprogramming] Niven's constant To: [email protected] Date: Wednesday, July 6, 2016, 12:40 AM Hi, I don't have my PC at hand now, but what is the problem with your calculation? I had never heard of Nivens constant, but looking at Wikipedia there is a formula for it containing the zeta function. Why not use this formula? Zeta 2 3 4 are well known constants, and you can calculate or hardcode the values for other even arguments, https://en.m.wikipedia.org/wiki/Particular_values_of_Riemann_zeta_function You might need to calculate bernoulli numbers. For odd values I think you need to calculate from scratch, but since zeta 14 is almost 1 you don't need to do too many calculations, just from 5 to 15, say. This is just thinking out loud. But my main point is you only probably need the first 15 or so zeta values before convergence is 'good enough'. Regards, Jon On Wed, Jul 6, 2016 at 12:12 AM +0900, "mikel paternain" <[email protected]> wrote: ________________________________ De: mikel paternain Enviado: domingo, 3 de julio de 2016 12:07:37 Para: [email protected] Asunto: Niven's constant Hi, I need improve the code for calculating Niven's constant. Could you help me? A first test is (see "Computo ergo sum" section in JoJ): exp=:+/"1@=@q: max=:>./ a=.max@exp(2+1.1000) C=.(+/%#)a Thanks in advance, Mikel Paternain ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
