[VA] SRC#001  Spiky Integral

07152018, 09:03 PM
Post: #21




RE: [VA] SRC#001  Spiky Integral
.
Hi, Gerson: Gerson W. Barbosa Wrote:I haven't trodden any of the hard (and beautiful) paths you all have done. Instead, I took an unallowed (according to Valentin's rules) and dull shortcut. I never said any of that, Gerson, and I didn't stablish any rules whatsoever so you hardly could've taken "an unallowed { VA: by whom ? not me ... } and dull shortcut." The only thing I said was this:
Thanks again for your long and very interesting recent post detailing the process you followed to get your results, that's exactly the spirit and the way for interested people to learn interesting things but in the future please do not credit me with statements, whether direct or implied, which I didn't actually make. Regards. V. . 

07152018, 09:22 PM
Post: #22




RE: [VA] SRC#001  Spiky Integral
(07152018 07:53 PM)DavidM Wrote: Gerson's method above is also how I would do it if the ListExt library isn't available. Or equivalently, using only builtin commands: << 0 m NDUPN >LIST >> 

07152018, 09:55 PM
Post: #23




RE: [VA] SRC#001  Spiky Integral
(07152018 09:03 PM)Valentin Albillo Wrote: . Hello, Valentin, Sorry for the misunderstanding! No, that was not what has elicited my comment about my possible infringement to your rules, which basically are “not googling for solutions”. Since I did google for an OEIS sequence, which eventually led me to a solution, it later appeared to me that I had done something that might be considered “illegal” or “unallowed”. Thanks for your always interesting problems and solutions! Gerson. 

07152018, 10:16 PM
Post: #24




RE: [VA] SRC#001  Spiky Integral  
07162018, 12:31 AM
(This post was last modified: 07162018 12:53 AM by Gerson W. Barbosa.)
Post: #25




RE: [VA] SRC#001  Spiky Integral
(07152018 06:26 PM)Thomas Klemm Wrote: This RPL program calculates the coefficients: Thanks! We now can use it to compute the integrals much faster than I was able to before: %%HP: T(3)A(R)F(.); \<< DUPDUP DUPDUP * + 2 / 2 MOD IF NOT THEN \<< \<< \> m \<< 0 m NDUPN \>LIST \>> \>> \> n ZEROS \<< { 0 1 } 2 n FOR k DUP SIZE \> a s \<< k ZEROS EVAL a + a 1 k SUB 0 + REVLIST s 1  ZEROS EVAL + ADD a k 1 + s SUB k 2 * ZEROS EVAL + ADD \>> NEXT \>> \>> EVAL 1 GET 2 ROT 2  ^ / \pi * ELSE DROP2 0 END \>> (I've only replaced { 1 m } 0 CON OBJ→ DROP m →LIST with John Keith's contribution elsewhere in this thread). 20 > '1909/65536*π', after 15.1 seconds (previously 118.4 seconds) And here is a list containing the results of the integrals, starting from n = 1 up to n = 71: { 0 0 '1/2*π' '1/4*π' 0 0 '1/8*π' '7/64*π' 0 0 '35/512*π' '31/512*π' 0 0 '361/8192*π' '657/16384*π' 0 0 '2055/65536*π' '1909/65536*π' 0 0 '24955/1048576*π' '46923/2097152*π' 0 0 '316301/16777216*π' '299973/16777216*π' 0 0 '4136805/268435456*π' '15796439/1073741824*π' 0 0 '13853361/1073741824*π' '26585247/2147483648*π' 0 0 '756388295/68719476736*π' '182188585/17179869184*π' 0 0 '20965992017/2199023255552*π' '20268008015/2199023255552*π' 0 0 '294245741167/35184372088832*π' '570497115729/70368744177664*π' 0 0 '4173319332859/562949953421312*π' '4055330794367/562949953421312*π' 0 0 '59723919552183/9007199254740992*π' '58153763705741/9007199254740992*π' 0 0 '430665931945033/72057594037927936*π' '840170667413757/144115188075855872*π' 0 0 '12505857230438737/2305843009213693952*π' '12217503312833669/2305843009213693952*π' 0 0 '182650875111521033/36893488147419103232*π' '44670833701814021/9223372036854775808*π' 0 0 '335205518724079925/73786976294838206464*π' } Just a few minutes on the emulator. Thanks again for providing both the program and an explanation why this works! Gerson. PS: Here is the result for n = 100, in case someone wants to check it :) '432756001487181254158446581/158456325028528675187087900672*π' (399 seconds, on the emulator) 

07162018, 03:27 AM
Post: #26




RE: [VA] SRC#001  Spiky Integral
(07152018 09:22 PM)John Keith Wrote: Or equivalently, using only builtin commands:That's what I was originally looking for but it appears to be missing on the HP48. Thanks for all your suggestions. (07162018 12:31 AM)Gerson W. Barbosa Wrote: We now can use it to compute the integrals much faster than I was able to before. If you want to create a list of all values you better do that within the FORloop: Code: « Kind regards Thomas 

07162018, 01:18 PM
Post: #27




RE: [VA] SRC#001  Spiky Integral
(07162018 03:27 AM)Thomas Klemm Wrote:(07152018 09:22 PM)John Keith Wrote: Or equivalently, using only builtin commands:That's what I was originally looking for but it appears to be missing on the HP48. My mistake, I was assuming you were using the HP 50g. 

07182018, 12:36 AM
Post: #28




RE: [VA] SRC#001  Spiky Integral
.
Hi, Gerson: Gerson W. Barbosa Wrote:PS: Here is the result for n = 100, in case someone wants to check it :) Checked Ok. I first get this: I(100) = 0.00857992304708725528996213203349184868757525160220371770625362449814379411183518043218667675290830524763068668295630334627497296 which then gets recognized as: 432756001487181254158446581 / 158456325028528675187087900672 * Pi in perfect agreement with your result. In return, will you please confirm my result for N=200, namely: I(200) = 195115902556687929766460554749767560813889646699192346811 / 200867255532373784442745261542645325315275374222849104412672 * Pi Regards. V. . 

07182018, 03:07 AM
Post: #29




RE: [VA] SRC#001  Spiky Integral
(07182018 12:36 AM)Valentin Albillo Wrote: In return, will you please confirm my result for N=200, namely: I am sorry, but I fear the HP 50g is not up to this task (at least with the current algorithm ). I(100) has required a clean emulated 50g and almost seven minutes, so I won’t even try. Anyway, the first seven significant digits of my approximate result on Free42, 3.05164078210624e3, agree with your numerical result for I(200). Not meaning to abuse your good will, would you please check how many significant digits I get right for N=20000? 3.06979593309e6 Thank you very much! Gerson. 

07182018, 06:17 AM
Post: #30




RE: [VA] SRC#001  Spiky Integral
For those who want NDUPN on their 48:
straightforward Code: \<< \> N \<< 0. N START DUP NEXT DROP2 N \>> \>> without local variables Code: \<< 0. OVER START OVER SWAP NEXT ROT ROT DROP2 \>> fast for large N Code: \<< Take your pick ;) Cheers, Werner 

07182018, 05:32 PM
(This post was last modified: 07182018 05:35 PM by ijabbott.)
Post: #31




RE: [VA] SRC#001  Spiky Integral
Is there a neat formula for just the constant term when converting the product to a sum? I tried Thomas Klemm's `spiky` function for N=20001 and it didn't like it very much! (Probably because it stores the coefficients for all the sequences up to N, eating lots of memory in the process.)
— Ian Abbott 

07182018, 11:39 PM
Post: #32




RE: [VA] SRC#001  Spiky Integral
.
Hi, Gerson: (07182018 03:07 AM)Gerson W. Barbosa Wrote: Not meaning to abuse your good will, would you please check how many significant digits I get right for N=20000? No, sorry, I can't verify your alleged numeric result for N=20,000, because:
Needless to say, computing I(20000) this way would require a truly humongous amount of time, certainly it would for my POPS (Pretty Old Pretty Slow) system, and regrettably I can't allocate that much running time to this task. Anyway, on a more feasible scale and in case it still might be useful to you, this is what a sufficiently accurate algorithm should return for N=1,000: I(1000) = 0.0002742581536 (all digits shown are correct) Regards. V. . 

07192018, 12:27 AM
Post: #33




RE: [VA] SRC#001  Spiky Integral
(07182018 05:32 PM)ijabbott Wrote: Is there a neat formula for just the constant term when converting the product to a sum? From A058377: Quote:FORMULA Thus I doubt there is a "neat formula". However here's a table of n, a(n) for n = 1..3342 From this we can calculate the value of I(1000) in accordance with Valentin's result: Code: >>> pi * 233854293495526890065238464248235751245240100131956760799187425375936277246905222228775610040031305849064670042246460679001271435777190838367631750031993315261824635725695903178339381850276654476951936755781933577607867005498666404483087885857942123647648138674089597298681941927332215904466338708 / 2**998 Since I had trouble with the Python program with bigger numbers I implemented it in Clojure: Code: (defn zeros [k] It's fast for n=100, takes a couple of seconds for n=200 and multiple minutes for n=1000. I've tried it for n=2000 and after a while that felt like eternity the correct result was given. From these measurements I would assume that it would take days or even weeks to calculate the value for n=20,000. Thus I refrained from trying. Kind regards Thomas 

07192018, 01:22 AM
Post: #34




RE: [VA] SRC#001  Spiky Integral
(07192018 12:27 AM)Thomas Klemm Wrote:(07182018 05:32 PM)ijabbott Wrote: Is there a neat formula for just the constant term when converting the product to a sum? There is an asymptotic formula, but it doesn't help much: https://cs.uwaterloo.ca/journals/JIS/VOL...ivan8.html (((n^2+n)/2+1) mod 2)*sqrt(6/pi)*2^(n1)/(n*sqrt(n)) n = 3 > 1.06385 (1) n = 4 > 1.38198 (1) n = 8 > 7.81764 (7) n = 11 > 38.7893 (35) n = 27 > 661051 (632602) n = 1000 > 2.34135e296 (2.3385429e296) A small correction might help a bit: (((n^2+n)/2+1) mod 2)* sqrt(6/pi)*2^(n1)*(16/(5*n)+21/(20*n^2)1/(8*n^3)+3/n^4)/(n*sqrt(n)) n = 3 >0.7969 (1) n = 4 > 1.07157 (1) n = 8 > 6.77707 (7) n = 11 > 34.8986 (35) n = 27 > 632623 (632602) n = 1000 > n = 1000 > 2.33854293231e296 (2.33854293496e296) Regards, Gerson. 

07192018, 01:31 AM
Post: #35




RE: [VA] SRC#001  Spiky Integral
(07182018 11:39 PM)Valentin Albillo Wrote: Needless to say, computing I(20000) this way would require a truly humongous amount of time, certainly it would for my POPS (Pretty Old Pretty Slow) system, and regrettably I can't allocate that much running time to this task. It has been useful! On Free42 I get 0.000274258153298, which is good to 9 significant digits. Since that's based on an asymptotic formula, the result for N=20,000 should be even more accurate. Thanks again, Gerson. 

07192018, 09:07 PM
Post: #36




RE: [VA] SRC#001  Spiky Integral
.
Hi, Gerson: (07192018 01:31 AM)Gerson W. Barbosa Wrote:(07182018 11:39 PM)Valentin Albillo Wrote: Anyway, on a more feasible scale and in case it still might be useful to you, this is what a sufficiently accurate algorithm should return for N=1,000: Seems likely. As I said in some previous post, I don't use any theoretical way to compute the numerators S(n) which, when divided by the respective powers of 2 (and times Pi), directly give the value of the integral. I simply compute the integral itself numerically using a quadrature algorithm, which for the case N=1,000 goes as follows (9 iterations): 0 0.000319732675251708806710904860429290813827540315 1 0.000293664662073707698824357734758191408742496989 2 0.000272327835887749444493527986927062286291843477 3 0.000274034638824020953737856637483929708504259670 4 0.000274259426644445675375095867663504951685809235 5 0.000274258154414552357469096400753858778167613322 6 0.000274258153608375557632839196888728683604929730 7 0.000274258153608378926807741119028301299681276600 8 0.000274258153608378926807734432669808007752908528 9 0.000274258153608378926807734432669808007979394750 So we get I(1000) = 0.000274258153608378926807734432669808007979394750 (all 48 decimal digits shown are correct) Applying a multilevel extrapolation scheme to those iterations quickly gives in excess of 100 correct decimal digits. These quadratureprovided results are useful to check that the S(n)provided ones are correct. Regards. V. 

07202018, 02:30 PM
(This post was last modified: 07202018 09:41 PM by Albert Chan.)
Post: #37




RE: [VA] SRC#001  Spiky Integral
(07102018 10:10 PM)Valentin Albillo Wrote: Hi all, welcome to my SRC#001  Spiky Integral: It seems the problem is easier to solve with symmetry. let F = cos(x) cos(2x) cos(3x) ... cos(Nx) I(N) = \(\int_{0}^{2 \pi} F dx \) = 2 \(\int_{0}^{\pi} F dx \) = 2 \(\int_{0}^{\pi/2} F dx \) + 2 \(\int_{\pi/2}^{\pi} F dx \) If F is symmetric around x = Pi/2, the two terms are same sized and same sign. A simple test is when x=Pi, F=1, which imply even number of odd values between 1 to N: (since cos(odd Pi) = 1, even numbers of such terms restored the symmetry) > If N mod 4 = 0 or 3, I(N) = 4 \(\int_{0}^{\pi/2} F dx \) Since F is maximized (= 1.0) when x=0, above should always be positive. As N increases, F spike is "thinner", thus smaller I(N), but still above 0 If odd number of odd values between 1 and N, symmetry is flipped. The two terms still have same size, but opposite sign, cancelling each other. > If N mod 4 = 1 or 2, I(N) = 0 So, for interval [1, 10] and nonzero I(N), N = 3,4,7,8 

07202018, 08:06 PM
(This post was last modified: 07212018 03:41 AM by Albert Chan.)
Post: #38




RE: [VA] SRC#001  Spiky Integral
(07182018 03:07 AM)Gerson W. Barbosa Wrote: would you please check how many significant digits I get right for N=20000? Let F= cos(x) cos(2x) cos(3x) ... cos(N x) Since N mod 4 = 20000 mod 4 = 0, I(N) = \(\int_{0}^{2 \pi} F dx \) = 4 \(\int_{0}^{\pi/2} F dx \) For big N, integral is dominated mostly by the area of spike: I(N) ~ 4 \(\int_{0}^{\pi /(2N)} F dx \) I did the integral in Python (plain float): I(20000) ~ 4 * 7.67448983276e07 = 3.0697959331e06 Both values agreed each other, to 11 digits. Comment: it is not necessary to sum the full spike area. For x = Pi / (20N) (one tenth of spike base), F = 1.547e36, which contribution little to the sum. With this tighter base, I(20000) still converge to the same value, but only take 16 sec (instead of 145 sec) BTW, my computer is 20+ years old Dell P3, modern computer may only take few seconds. 

07212018, 04:06 AM
Post: #39




RE: [VA] SRC#001  Spiky Integral
(07202018 08:06 PM)Albert Chan Wrote: For big N, integral is dominated mostly by the area of spike: Thanks both for the I(20000) result and for another interesting approach. These results might be useful to improve the correction terms of the approximation formula. 

07212018, 04:00 PM
Post: #40




RE: [VA] SRC#001  Spiky Integral
(07212018 04:06 AM)Gerson W. Barbosa Wrote: Thanks both for the I(20000) result and for another interesting approach. These results might be useful to improve the correction terms of the approximation formula. Hi Gerson, I am a bit late to the game, and just realized your number is also an approximation. So, both numbers matching are nice, but probably not good enough. Instead, I tried I(1000) ~ 4 \(\int_{0}^{\pi /2000} F dx \) = 0.00027425815360837926 From Valentin Albillo last post, my estimate match correct values, to 15 digits Comparing values in binary form, accuracy is even better, mine is just 6 ULP over ! Thanks Valentin. I am amazed at how your integration function work. You mentioned the problems of it computing I(20000). Shrinking the integral range 800,000 times may help: (0 to Pi/400000) instead of (0, 2 Pi) 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)