On Sat, Jan 31, 2009 at 11:40 AM, William Stein <[email protected]> wrote:
> On Sat, Jan 31, 2009 at 11:26 AM, mabshoff
> <[email protected]> wrote:
>>
>>
>>
>> On Jan 31, 11:17 am, William Stein <[email protected]> wrote:
>>> On Sat, Jan 31, 2009 at 11:15 AM, Rolandb <[email protected]> wrote:
>>>
>>> > Hi,
>>>
>>> > I received first a MemoryError, and later on Sage reported:
>>>
>>> That means that Sage ran out of memory.  It's not a bug -- it's just a
>>> limitation of your computer.
>>
>> srange() builds a list in memory, so using an iterator will make that
>> code use less memory.
>
> Good idea.  His an iterator version of the code:
>
> uitkomst1=[]
> uitkomst2=[]
> eind=int((10^9+2)/(2*sqrt(3)))
> print eind
> for y in (1..eind):
>  test1=is_square(3*y^2+1,True)
>  test2=is_square(48*y^2+1,True)
>  if test1[0] and test1[1]%3==2: uitkomst1.append((y,(2*test1[1]-1)/3))
>  if test2[0] and test2[1]%3==1: uitkomst2.append((y,(2*test2[1]+1)/3))
> print uitkomst1
> een=sum(3*x-1 for (y,x) in uitkomst1 if 3*x-1<10^9)
> print uitkomst2
> twee=sum(3*x+1 for (y,x) in uitkomst2 if 3*x+1<10^9)
> print een+twee
>
> It prints
>
> 288675135
> ... and I got bored waiting for the rest...

OK, I got really tired of waiting, and just wrote a cython version of
the above, which does it for 10^9 in 25 seconds.  It's over 200 times
faster than your code.  See attached worksheet or paste this into a
notebook cell and press shift enter, then do f(10^9):


%cython
from sage.all import is_square
cdef extern from "math.h":
    long double sqrtl(long double)

def f(n):
    uitkomst1=[]
    uitkomst2=[]
    cdef long long eind=int((n+2)/(2*sqrt(3)))
    cdef long long y, t
    print eind
    for y in range(1,eind):
        t = <long long>sqrtl(<long long> (3*y*y + 1))
        if t * t == 3*y*y + 1:
            uitkomst1.append((y, (2*t-1)/3))
        t = <long long>sqrtl(<long long> (48*y*y + 1))
        if t * t == 48*y*y + 1:
            uitkomst2.append((y, (2*t+1)/3))
    print uitkomst1
    een=sum([3*x-1 for (y,x) in uitkomst1 if 3*x-1<10^9])
    print uitkomst2
    twee=sum([3*x+1 for (y,x) in uitkomst2 if 3*x+1<10^9])
    print een+twee

I don't guarantee there isn't some funny overflow or wrongness with
the above.  All I claim is that it is faster than you'll get use any
interprter.

William

--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to 
[email protected]
For more options, visit this group at 
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---

Attachment: sage_support_Rolandb.sws
Description: Binary data

Reply via email to