Sums of all Squares

[bsquared]

Originally Posted by alexhiggins732

n Pmax Sum(P^2) from 2 to Pmax; a multiple of 10^n
1 907 37464550
2 977 46403000
3 977 46403000
You have duplicated 2 and 3 on your site.

That's not an error... see post 8.

Originally Posted by alexhiggins732

I also noted that the OEIS sequence didn't have example code. Perhaps you could submit the expression you have entered into YAFU.

I could submit the PARI code others have used, for instance CRGreathouse:

 
s=Mod(0,10^4);forprime(p=2,1e9,s+=p^2;if(!s,return(p)))

or axn, in post 37:

p=2;s=0;ten=10;for(i=2,10000000,s=s+p;if(Mod(s,ten)==0, ten=ten*10;print(i,":",p, ":",s)); p=nextprime(p+1))

The YAFU code is a little too unwieldy for OEIS... it is a hundred lines or so of C code mixed in with some custom assembly and leveraging thousands of lines of more C code implementing the sieve of eratosthenes. I have no problem making it public, but wasn't aware of any interest in it.

[alexhiggins732]

Originally Posted by bsquared

That's not an error... see post 8.


I could submit the PARI code others have used, for instance CRGreathouse:

 
s=Mod(0,10^4);forprime(p=2,1e9,s+=p^2;if(!s,return(p)))

or axn, in post 37:

p=2;s=0;ten=10;for(i=2,10000000,s=s+p;if(Mod(s,ten)==0, ten=ten*10;print(i,":",p, ":",s)); p=nextprime(p+1))

The YAFU code is a little too unwieldy for OEIS... it is a hundred lines or so of C code mixed in with some custom assembly and leveraging thousands of lines of more C code implementing the sieve of eratosthenes. I have no problem making it public, but wasn't aware of any interest in it.

The more examples the merrier. Personally, I feel code examples help me to understand the sequences better, as I think can look at code an understand it easier then an explanation using a mathematical explanation.

The best example may be I frankly could not understand the first mathematical notation I read to explain the Lucas-Lehmer test, but looking at the code it understood the notation that was used immediately.

As far as the YAFU code, I was under the impression that you a simple PARI like expression was entered into YAFU.

I am also interested in how all primes up to 10 trillion have been calculated, or have I read it wrong and that is the sum of the primes is over 10 Trillion.

[bsquared]

Originally Posted by alexhiggins732

The more examples the merrier. Personally, I feel code examples help me to understand the sequences better, as I think can look at code an understand it easier then an explanation using a mathematical explanation.

The best example may be I frankly could not understand the first mathematical notation I read to explain the Lucas-Lehmer test, but looking at the code it understood the notation that was used immediately.

As far as the YAFU code, I was under the impression that you a simple PARI like expression was entered into YAFU.

YAFU is nowhere near as powerful an environment as PARI. The only reason I used it at all was because it has a rather fast prime generation routine. But I didn't "use" it by constructing some combination of built-in functions. It probably wasn't clear, but I "used" it by writing new code which leveraged functions elsewhere in the code base.

The way it's written, I doubt the YAFU code would make anything clearer to you .

Originally Posted by alexhiggins732

I am also interested in how all primes up to 10 trillion have been calculated, or have I read it wrong and that is the sum of the primes is over 10 Trillion.

Going on 100 trillion now, actually .

The sum is a number with 26 digits, and the square sum is a number with 39 digits. The primes are calculated using the sieve of eratosthenes. That code is public domain... just download YAFU-1.18 and browse around through the soe.c file. Although it's not well documented .

I compute the primes in batches of 1 billion integers. At a height of 10 trillion, that is 32 million primes or so per batch. It's then pretty straightfoward to sum them and check each for divisiblity by a power of 10 (sweeping many optimizations under the rug, for clarity). On my (admittedly fast) computer, each batch takes about 3 seconds for prime computation, and .2 seconds for the sum/check. This is the only code that's not public yet. Although yafu-1.19 will probably retain it when it comes out.

[alexhiggins732]

Originally Posted by bsquared

YAFU is nowhere near as powerful an environment as PARI. The only reasons I used it at all was because it has a rather fast prime generation routine. But I didn't "use" it by constructing some combination of built-in functions. It probably wasn't clear, but I "used" it by writing new code which leveraged functions elsewhere in the code base.



The way it's written, I doubt the YAFU code would make anything clearer to you .

I do not seek clarification on this, I understand this. I was making a generalized statement and further it may not be clear immediately reading the description of the sequence. However, a code example helps with that clarification.

Here's .Net Code If you wish to post as well.

VB

dim ten as bigint=10:dim sum as bigint=2:for each prime in primesbelow(10000000):if ten.mod(sum)=0 then : ten*=10 : console.writeline("{0}:{1}:{2}",i,p,s): end if: sum+=prime: next

C#

 
bigint ten = 10;bigint sum = 2; foreach (var prime in  primesbelow(10000000)) { if (ten.mod(sum) == 0) {ten *= 10; console.writeline("{0}:{1}:{2}", i, p, s); } sum += prime; }}

[alexhiggins732]

Originally Posted by bsquared

YAFU is nowhere near as powerful an environment as PARI. The only reason I used it at all was because it has a rather fast prime generation routine. But I didn't "use" it by constructing some combination of built-in functions. It probably wasn't clear, but I "used" it by writing new code which leveraged functions elsewhere in the code base.

The way it's written, I doubt the YAFU code would make anything clearer to you .



Going on 100 trillion now, actually .

The sum is a number with 26 digits, and the square sum is a number with 39 digits. The primes are calculated using the sieve of eratosthenes. That code is public domain... just download YAFU-1.18 and browse around through the soe.c file. Although it's not well documented .

I compute the primes in batches of 1 billion integers. At a height of 10 trillion, that is 32 million primes or so per batch. It's then pretty straightfoward to sum them and check each for divisiblity by a power of 10 (sweeping many optimizations under the rug, for clarity). On my (admittedly fast) computer, each batch takes about 3 seconds for prime computation, and .2 seconds for the sum/check. This is the only code that's not public yet. Although yafu-1.19 will probably retain it when it comes out.

Must be very fast. I use a "similar sliding window" technique if you wish to call it that. However, after each window I need to save the primes to a .bin file in the form of 6k+=1, (eg, skip 2 and 3 and 5,7,11,13,17,19,23,25 =11110000, for each bit and after each iteration start at the beginning of the binary file and calculate the start index of each prime upto and sieve the current window, saving results to binary and and continuing.

However, for all primes below 2^32 the file size is 178,956,968 bytes (170 MB) and it doubles with each power of two which led me to ask how how you could sieve up to 2^44 which would require a huge amount of memory even when saved in the most compact possible form (1 bit for each 6K+=1)

Which brings up another point. If the distribution of primes are indeed random, how then is it possible to compress the aforementioned file and receive a compression ratio of 60% (eg 170mb compresses to 102MB)?

[CRGreathouse]

If you use my code, be sure to change "4" to "n" (and perhaps prefix "a(n)=").

Originally Posted by bsquared

The YAFU code is a little too unwieldy for OEIS... it is a hundred lines or so of C code mixed in with some custom assembly and leveraging thousands of lines of more C code implementing the sieve of eratosthenes. I have no problem making it public, but wasn't aware of any interest in it.

I would like to see it, actually. There are plenty of such submissions to the OEIS -- short programs in the %o (or %p or %t, of course) lines, long programs as links to source. But since you have a page set up, why not put it there?

[CRGreathouse]

Originally Posted by alexhiggins732

However, for all primes below 2^32 the file size is 178,956,968 bytes (170 MB) and it doubles with each power of two which led me to ask how how you could sieve up to 2^44 which would require a huge amount of memory even when saved in the most compact possible form (1 bit for each 6K+=1)

Why save the primes? Generate a batch, do what you need to do, then discard them. The only primes you need to keep from iteration to iteration are the sieving primes. If you want to test up to 100 trillion, that's just the first 664,579 (2.5 MB using 32-bit ints). A bit array of a billion integers takes a billion bits, or 119 MB.

[bsquared]

Originally Posted by CRGreathouse

Why save the primes? Generate a batch, do what you need to do, then discard them. The only primes you need to keep from iteration to iteration are the sieving primes. If you want to test up to 100 trillion, that's just the first 664,579 (2.5 MB using 32-bit ints). A bit array of a billion integers takes a billion bits, or 119 MB.

Exactly. Although I don't even keep the sieving primes because that's overhead I don't want to manage and its only a hundredth of a second to compute them. Using -v -v prints some more info in yafu-1.18. For example, here is a yafu command to count the primes in a range from 1e13 to 1.0001e13:

yafu "primes(10000000000000,10001000000000)" -v -v

which yeilds this output:

elapsed time for seed primes = 0.0144
sieving range 9999999999840 to 10001006632800
using 227655 primes, max prime = 3162436
using 8 residue classes
lines have 4194304 bytes and 33554432 flags
lines broken into 128 blocks of size 32768
blocks contain 262144 flags and cover 7864320 primes
bucket sieving 194320 primes > 393216
allocating space for 62341 hits per bucket
using 45026396 bytes for sieving storage
elapsed time = 1.0610
ans = 33405006

Counting this range only takes a second. Even with bucket sieving, all told it only uses 45 MB of RAM, because the sieve doesn't store a bit for every number, just the ones not divisible by 2, 3, or 5 (i.e. I line sieve the residue classes mod 30).

Once I know the locations of the primes in the bit arrays, I compute them, but only the 33 odd million of them, of course, not a billion. Stored as 64 bit integers, this takes about 250 MB of RAM and another 2 seconds or so.

[bsquared]

Originally Posted by CRGreathouse

I would like to see it, actually. There are plenty of such submissions to the OEIS -- short programs in the %o (or %p or %t, of course) lines, long programs as links to source. But since you have a page set up, why not put it there?

Thanks for your interest! I did have fun with the code for this puzzle, so it's neat to be able to share it. I haven't put it on my site yet because it's not useable as a standalone program due to the calls to internal yafu functions. I don't know if its worth the effort to repackage it as a standalone program, but it will be available within yafu once I release a new version.

To save on the forum database, its attached rather than posted, but keep in mind that this won't compile as is.

- ben.

[bsquared]

The routine to do prime cubes is very similar, except for the assembly to do the cube/sum:

ASM_G (
"movq %1, %%rcx \n\t"               /* store prime */
"mulq %%rcx \n\t"                   /* square it */
"movq %%rax, %%r8 \n\t"           /* save p^2 lo (a) */
"movq %%rdx, %%r9 \n\t"           /* save p^2 hi (d) */
"mulq %%rcx \n\t"                   /* p * a */
"movq %%rax, %%r10 \n\t"         /* save p*a lo (apa) */
"movq %%rdx, %%r11 \n\t"         /* save p*a hi (apd) */
"movq %%r9, %%rax \n\t"           /* p * d */
"mulq %%rcx \n\t"                  /* lo part in rax (dpa), hi in rdx (dpd) */
"addq %%r10, (%%rbx) \n\t"       /* sum0 = sum0 + apa */
"adcq %%rax, 8(%%rbx) \n\t"      /* sum1 = sum1 + dpa + carry */
"adcq %%rdx, 16(%%rbx) \n\t"    /* sum2 = sum2 + dpd + carry */
"addq %%r11, 8(%%rbx) \n\t"      /* sum1 = sum1 + apd */
"adcq $0, 16(%%rbx) \n\t"          /* sum2 = sum2 + carry */
: 
: "b"(sum->val), "a"(PRIMES[j])
: "rcx", "rdx", "r8", "r9", "r10", "r11", "memory", "cc"); 

[alexhiggins732]

Originally Posted by bsquared

YAFU

Just download YAFU-1.18 and browse around through the soe.c file. Although it's not well documented .

I compute the primes in batches of 1 billion integers. At a height of 10 trillion, that is 32 million primes or so per batch. It's then pretty straightfoward to sum them and check each for divisiblity by a power of 10 (sweeping many optimizations under the rug, for clarity). On my (admittedly fast) computer, each batch takes about 3 seconds for prime computation, and .2 seconds for the sum/check. This is the only code that's not public yet. Although yafu-1.19 will probably retain it when it comes out.

Downloaded the source and fired up the Visual Studio Debugger after a bit of .... and still a bit hard to see exactly what is going on because pointers to the arrays don't allow you to see the contents of the array on the memory address. While c/c++ is second language I can't say I fully understand pointer arithmetic entirely, but I am sure if I could figure it out.

So being the crank I am, I did what any self-respecting crank would do, I rewrote it in VB.NET, but have a few problems perhaps you can shine some light on.

I tried, for the most part, to use the exact vb equivalent of the c code
1) The multi-threading is not functional, yet. I understand the logic that you are implementing so I just may rewrite that and deviate from adhering strictly to the original code. I couldn't really figure out what thread_soedata_t->thread_id was being used for. I am not sure what the vb equivalent of CreateThread would be. In any case, I can rewrite the threading. Each thread is assigned a command, waits until it is signal to run then notifies the master thread when done, rinse lather repeat.

2) When you declare a local variable of a structure member a value copy is created on the heap and the original structure is not updated. If this where declared a class, it would work like the c implementation but I didn't what to do that. So I had to rewrite as follows:

Public Sub primes_from_lineflags(ByRef t As thread_soedata_t)
         ' this is no good. it create value copy on the heap so the original
         ' thread_soedata_t members do not get updated.
         'Dim ddata As soe_dynamicdata_t = t.ddata
        'Dim sdata As soe_staticdata_t = t.sdata
        ' work around is to pass a pointer to the structure and access the
        ' the members using the pointer like on the next line.
        Dim line() As Byte = t.ddata.line

Problem is, I think that there is somewhere else in the code the value copy of the structure is being accessed and I can't figure out where. Stepping through the yafu source code and the vb source code everything matches up except, the number of primes that are returned. Only thing I can think of is due to pointer arithmetic some values the pointer arrays are different than the vb version.

However, not being able to see that values are in those arrays in the debugger and being a crank that is trying to learn what is going on from the code I don't know if the values in the vb arrays are correct.

Maybe my misinterpretation of how pointers work is the problem.

For example, interpreted flagblock += BLOCKSIZE; as increment the index of the flagblock + blocksize. Since this is illegal in vb, to accommodate I created a new variable fob and instead of flagblock + BLOCKSIZE, its fob+BLOCKSIZE and then when accessing flagblock its flagblock[x+fob]. Is this wrong?


Perhaps the simplest solution is to attach the code I have and hopefully someone with MSVS can debug it and give me a hint whats wrong. Its just a little frustrating after rewriting 2000 lines of code....

Oh well, maybe I'll write it in C# so I can use the pointer arithmetic and just import the DLL into my vb projects