Re: Re: Re: Re: Re: Timing
- To: mathgroup at smc.vnet.net
- Subject: [mg88087] Re: [mg88069] Re: [mg88017] Re: [mg88011] Re: [mg87971] Re: Timing
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 24 Apr 2008 05:56:57 -0400 (EDT)
- References: <480A0EE6.3090109@csl.pl> <200804211038.GAA28042@smc.vnet.net> <200804211840.OAA09974@smc.vnet.net> <200804221026.GAA28529@smc.vnet.net> <200804230810.EAA29058@smc.vnet.net>
Artur wrote:
> Many thanks for Andrzej Kozlowski for this procedure!
> I was execute these on my computer and received follwing Timing
> Table[f[n], {n, 2, 20}] // Timing = 0.071
> Table[f[n], {n, 2, 21}] // Timing = 0.089
> Table[f[n], {n, 2, 22}] // Timing = 525.756
> but increasing limit of n to 23 was catastrophe on my computer (probably
> RAM memory is limit)
>
> Could help me sombody with big memory size computer check these
> procedure for n=23 (next true should occured for n=31)
> Or sombody can explained this catastrophe in relation hardware-software
> relation or estimated how increased time if increased n (is this time
> polynomial time or not)
>
> Best wishes
> Artur
Apparently you wish to replicate
http://www.research.att.com/~njas/sequences/?q=3%2C5%2C7%2C13%2C17%2C19%2C31%2C61%2C89&language=english&go=Search
First, the way you are going about this is hopeless. You are generating
huge integers, on the order of many millions of bits, just to see if
they are divisible by much smaller integers on the order of perhaps
dozens of bits. This "catastrophe" is entirely self inflicted; you
simply MUST avoid the huge numbers in the first place. Functions like
Mod can be real helpful, in such circumstances.
So you can tackle it as follows. Treat Sqrt[2] as a variable, that is,
use polynomials in x, reduced by x^2-2. Also you want to do reduction
modulo 2^n-1, since all you care about is whether the remainder is zero.
This can be interleaved with powering, for efficiency (avoids huge
numbers/expressions).
Back before version 6 there was an add-on package PolynomialPowerMod
that handled such things, but (I think) was restricted to prime number
moduli. Anyway, I got a bit frustrated trying to figure out how to do
this efficiently using built-in functions in version 6. This is an
admission that either the documentation was lacking, or I was too tired.
Also, I ran afoul of a difficult-to-track crash bug in our development
kernel, and I need to look into that. So I just wrote a bit of code to
do what we need here.
ppMod[poly_, n_, {p2_,m_}] := Module[
{bits=Reverse[IntegerDigits[n,2]], pj=poly, oldpj, res=1},
If [bits[[1]]==1, res = pj];
Do[
pj = PolynomialMod[pj^2, {p2,m}];
If [bits[[j]]==1, res = PolynomialMod[res*pj,{p2,m}]],
{j, 2, Length[bits]}];
res
]
Timing[Map[First,Select[Table[{n, PolynomialMod[
(ppMod[3+2*x,2^(n-1)-1,{x^2-2,2^n-1}] -
ppMod[3-2*x,2^(n-1)-1,{x^2-2,2^n-1}])*(x/8),
{x^2-2,2^n-1}]}, {n,1,1000}], #[[2]]===0&]]]
Out[10]= {158.978, {1, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607}}
Daniel Lichtblau
Wolfram Research
- References:
- Re: Timing
- From: Artur <grafix@csl.pl>
- Re: Re: Timing
- From: "W_Craig Carter" <ccarter@mit.edu>
- Re: Re: Re: Timing
- From: Andrzej Kozlowski <akoz@mimuw.edu.pl>
- Re: Re: Re: Re: Timing
- From: Artur <grafix@csl.pl>
- Re: Timing