some math behind it, for 1e5 S y, 1763199 is the maximum y for which 11 is a sufficient prime.
q: >: 1e5 -~ 1763200 + 7559 2 2 2 3 3 3 5 7 13 17 this is a d of 256, but boosts the S total by 4, and requires 13. For this difference of 4, the number of d 1e5 less than the max range must be higher than 256, and the next highest d within the 1e5 range must be 252 (4 less than 256) q: >: 1e5 -~ 1763200 + 50399 2 2 2 2 2 2 3 3 5 5 7 17 d >: 1e5 -~ 1763200 + 50399 252 (does not use 13) d 1e5 -~ 1763200 288 Its interesting that the range u can go beyond 4e8 without increasing 13, but from 4e8 to 1e9, it has to go up all the way to 31. I started this post with the statement that there was math... but maybe its luck. ----- Original Message ----- From: Mike Day <[email protected]> To: [email protected] Cc: Sent: Saturday, October 25, 2014 1:20 PM Subject: Re: [Jprogramming] An easy Euler Project one (485) Wow! 31 is a sufficient maximum for the full Problem 485; 29 isn't. (Trial and error, not mathematical analysis!!!) It takes about 27 seconds with 31 as maximum prime compared to the 45 seconds using all primes <: %: n . Cheers, Mike PS I played around with that "sliding window" article I cited earlier. Not surprisingly, it crawls like a snail with an effort to program it in J. Too loopy. Perhaps I could tacitise it, but not easy and hardly worth the effort. On 25/10/2014 17:51, 'Pascal Jasmin' via Programming wrote: > > combined our approaches with a simple modification to dsieve > > add 2 lines to top of function: > ( %: y) dsieve y > : > > change the while line to: > > while. p <: x do. > > for 1e5 S 4e8, you do not need to look at primes higher than 13, and so this > runs in 9 seconds > > timespacex '100000 ([: +/ >./\) 13 dsieve 40000000' > 9.74692 2.51659e9 > > Smaller ranges (y) increase the required prime relative to max u (x)... 1e4 S > 4e6 requires primes up to 23 > > 10000 ([: +/ >./\) 23 dsieve 4000000 > 932689064 > > There is probably a mathematical relationship but a simple conjecture would > be that if including the next prime in the sieve loop does not increase the S > total, then no other prime above that would either. So there could be a > built in loop short circuit. > > > ________________________________ > From: Mike Day <[email protected]> > To: [email protected] > Sent: Saturday, October 25, 2014 6:29 AM > Subject: Re: [Jprogramming] An easy Euler Project one (485) > > > FWIW, I've finally cracked the minute for the full problem > using the improved sieve below. It takes about 45 sec on > my laptop. > > 3 points: > 0) Problem 485 solvers' discussion reveals this discussion about > maxima in a sliding window: > http://richardhartersworld.com/cri/2001/slidingmin.html > > I doubt if it outperforms J's elegant method, eg > 10 >./\ v > but I might have a look. > > 1) Several solvers seemed to have used a sieve, others brute- > force, as I did originally. Some discussion of the windowing, > including (0) above. > > 2) I'm puzzled why this method of producing an arithmetic > progression of indices into an array > i =. <:+/\p#~<.n%p > is quicker than the more obvious: > i =. I.n$1{.~ - p > > Here's the improved sieve; I expect there will be gratuitous > line-throws again: > > dsieve =: 3 : 0 > > NB. SLIGHTLY FASTER WITH m THAN s IN PREVIOUS VN > > m =. 1#~n =. y NB. running product array > > d =. m NB. number of divisors array > > NB. every number i is a product of pij^qij for some distinct primes, pij > > NB. required di is the product of (qij+1) > > p =. 2 NB. seed first prime candidate > > while. p <: %:n do. > > np =. (<.n%p)<. <. p^npl =. <: >.p^.n > > pl =. <. p^i.npl+1 > > > NB. get q = powers of prime p in multiples of p > > NB. eg q for p=3 is 1 1 2 1 1 2 1 1 3 for n < 27^2 > > NB. cf 3 6 9 12 15 18 21 24 27 > > q =. (}: , >:@{:)@:($~ (p * #))^:(0>.npl-1)1 > > q =. q$~<.n%p NB. extend q for all n(%p) > > NB. i =. I.n$1{.~ - p NB. indices of all multiples of p > > NB. next line, alternative for getting index to amend s and d, > > NB. is FASTER! > > i =. <:+/\p#~<.n%p > > > m =. (<. (i{m) * q{ pl) i } m NB. embed p^q for all multiples of p > > d =. ( (i{d) * q + 1) i } d NB. adjust d with contributions from p^q > > > NB. first three sets of m & d > > NB. +-+------------------------------------------------------------+ > > NB. |2|1 2 1 4 1 2 1 8 1 2 1 4 1 2 1 16 1 2 1 4 1 2 1 8 1 2 1 4 1 2| > > NB. | |1 2 1 3 1 2 1 4 1 2 1 3 1 2 1 5 1 2 1 3 1 2 1 4 1 2 1 3 1 2| > > NB. +-+------------------------------------------------------------+ > > NB. +-+----------------------------------------------------------------+ > > NB. |3|1 2 3 4 1 6 1 8 9 2 1 12 1 2 3 16 1 18 1 4 3 2 1 24 1 2 27 4 1 6| > > NB. | |1 2 2 3 1 4 1 4 3 2 1 6 1 2 2 5 1 6 1 3 2 2 1 8 1 2 4 3 1 4| > > NB. +-+----------------------------------------------------------------+ > > NB. > +-+---------------------------------------------------------------------+ > > NB. |5|1 2 3 4 5 6 1 8 9 10 1 12 1 2 15 16 1 18 1 20 3 2 1 24 25 2 27 4 > 1 30| > > NB. | |1 2 2 3 2 4 1 4 3 4 1 6 1 2 4 5 1 6 1 6 2 2 1 8 3 2 4 3 1 8| > > NB. > +-+---------------------------------------------------------------------+ > > > NB. smoutput p;m,:d > > NB. wd'msgs' > > > p =. 4 p: p NB. next prime > > NB. next line, alternative for next prime, is too slow > > NB. p =. ((1+p+1 i.~p&}.) ) m > > end. > > > NB. final m is 1+i.n except for primes >: %: n > > NB. double d for every such prime > > d * >:(<>:@i.@# ) m > > ) > > > Perhaps I'll stop here, who knows! > > Mike > > On 23/10/2014 16:06, Mike Day wrote: >> In case anyone's interested, here's a faster, if slightly >> more opaque, version of a sieve for number of divisors >> for all positive integers up to a limit. >> >> It's better commented though. >> >> Sorry for the inevitable double line-throws. >> >> dsieve =: 3 : 0 >> >> s =. <.>:i.n=.y NB. original numbers to be "factored" >> >> d =. n#1 NB. number of divisors array >> >> NB. every number i is a product of pij^qij for some distinct primes, pij >> >> NB. required di is the product of (qij+1) >> >> p =. 2 NB. seed first prime candidate >> >> while. p <: %:n do. >> >> np =. (<.n%p)<.<. p^npl =. <: >.p^.n >> >> pl =. <. p^i.npl+1 >> >> >> NB. get q = powers of prime p in multiples of p >> >> NB. eg q for p=3 is 1 1 2 1 1 2 1 1 3 for n < 27^2 >> >> NB. cf 3 6 9 12 15 18 21 24 27 >> >> q =. (p&((}:,>:@{:)@(]$~(*#)))^:(0>.npl-1))1 >> >> >> q =. q$~<.n%p NB. extend q for all n(%p) >> >> i =. I.n$1{.~ - p NB. indices of all multiples of p >> >> NB. next line, alternative for getting index to amend s and d, >> >> NB. is slower >> >> NB. i =. p&*&.>:i.<.n%p >> >> s =. (<. (i{s) % q{ pl) i } s NB. remove p^q from multiples of p >> >> d =. ( (i{d) * q + 1) i } d NB. adjust d with contributions from p^q >> >> p =. 4 p: p NB. next prime >> >> NB. next line, alternative for next prime, is too slow >> >> NB. p =. (>:@i.&1 @: >&1 ) s >> >> >> NB. after p=2, s = 1 1 3 1 5 3 7 1 9 5 11 3 13 7 15 1 ... >> >> NB. after p=2, d = 1 2 1 3 1 2 1 4 1 2 1 3 1 2 1 5 ... >> >> >> NB. after p=3, s = 1 1 1 1 5 1 7 1 1 5 11 1 13 7 5 1 ... >> >> NB. after p=3, d = 1 2 2 3 1 4 1 4 3 2 1 6 1 2 2 5 ... >> >> >> NB. etc... >> >> >> NB. smoutput p;d,:s >> >> NB. wd'msgs' >> >> end. >> >> >> NB. final s is all ones except for primes >: %: n >> >> NB. double d for every such prime >> >> d * 1 2{~ s>1 NB. or d * >: s > 1 >> >> ) >> >> >> Ssieve =: (+/@:(>./\ dsieve )~)/ NB. solves problem 485 in about 5 >> minutes. >> >> Ssieve 1000 10 >> >> 17176 >> >> >> Mike >> >> On 22/10/2014 11:48, Mike Day wrote: >>> Having floated this boat, here's a sieve method for >>> obtaining all numbers of divisors for integers from >>> 1 to n: >>> >>> dsieve =: 3 : 0 >>> >>> s =. <.>:i.n=.y >>> >>> d =. n#1 >>> >>> p =. 2 >>> >>> while. p < %:n do. >>> >>> pp =. 1 >>> >>> np =. (<.n%p)<.<.p^npl =. <.p^.n >>> >>> q =. np#0 >>> >>> pl =. p^i.npl+1 >>> >>> NB. get q = powers of prime p >>> >>> NB. eg q for p=3 is 1 1 2 1 1 2 1 1 3 ... >>> >>> while. pp<n do. >>> >>> q =. q+np$1{.~-pp >>> >>> pp =. <.pp*p >>> >>> end. >>> >>> NB. q =. +/@:((+/) ($(1{.~-) )"0 ]) pl NB. slower than the explicit loop >>> >>> q =. q$~<.n%p >>> >>> i =. I.mask =. n$1{.~ - p >>> >>> s =. (<. (mask#s) % q{pl) i } s >>> >>> d =. ((mask#d) * q + 1) i } d >>> >>> p =. 4 p: p NB. or could search for next non-1 in }.s >>> >>> end. >>> >>> d * >:s>1 >>> >>> ) >>> >>> >>> timer'#dsieve 1000000' >>> >>> +-----+-------+ >>> >>> |1.446|1000000| >>> >>> +-----+-------+ >>> >>> >>> Pari GP does this in about 1.2 seconds, so a reasonable >>> >>> comparison. >>> >>> >>> Thanks, >>> >>> >>> Mike >>> >>> >>> On 22/10/2014 00:15, 'Pascal Jasmin' via Programming wrote: >>>> a complication with the approach is: >>>> >>>> */ 2 2 2 3 3 5 5 7 11 >>>> 138600 >>>> >>>> d */ 2 2 2 3 3 5 5 7 11 >>>> 144 >>>> >>>> and so for the ranges starting at 38600 to 63320 there is a greater >>>> maximum, and it also appears that for some multiples of 138600, the >>>> number of divisors exceeds (for part of the range) of multiples of >>>> 83160 >>>> >>>> a number that simplifies the process is: >>>> >>>> */ 2 2 2 3 5 7 11 >>>> 9240 >>>> >>>> as it allows needing to examine with d only multiples of it. >>>> >>>> >>>> >>>> ----- Original Message ----- >>>> From: Raul Miller <[email protected]> >>>> To: Programming forum <[email protected]> >>>> Cc: >>>> Sent: Tuesday, October 21, 2014 6:41 PM >>>> Subject: Re: [Jprogramming] An easy Euler Project one (485) >>>> >>>> See also http://oeis.org/A002182 >>>> >>>> Also Roger Hui's point about a sieve is probably a good one since this >>>> problem only needs you to consider prime factors less than 10000. >>>> Still, >>>> those 1229 prime factors multiplied by the 1e8 limit means we would >>>> need >>>> something on the order of several terabytes of memory for intermediate >>>> results if we were to store the entire sieve in the obvious way. So >>>> some >>>> extra work would be needed, and I'm not sure how that would pan out. --- This email is free from viruses and malware because avast! Antivirus protection is active. http://www.avast.com ----- No virus found in this message. Checked by AVG - www.avg.com Version: 2014.0.4765 / Virus Database: 4040/8449 - Release Date: 10/25/14 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
