mm-829
===
Subject: Re: Zero times a variable in numerical calculations
> 
> Chop[0.0 a] indeed works and returns nicely 0 however I still have some 
> problems when using this on Sparse Arrays.
> 
> I define:
> M = SparseArray[{{1, 1} -> a, {2, 2} -> 1.0, {3, 3} -> 1.0, {2, 3} -> b, \
{3, 
> 2} -> b, {_, _} -> 0}, {3, 3}];
> and
> T = SparseArray[{{1, 1} -> 1.0, {2, 2} -> -Sqrt[0.5], {3, 3} -> Sqrt[0.5], \

> {2, 3} -> Sqrt[0.5], {3, 2} -> Sqrt[0.5], {_, _} -> 0}, {3, 3}];
> I calculate T.M.T This should result in a sparse array with 3 nonzero 
> elements. (a,1-b and 1+b on the diagonal).
> However
> T.M.T=SparseArray[<5>, {3, 3}]
> (T.M.T)[[3,2]]=0. (0.707107 + 0.707107 b)
> and
> (T.M.T)[[2,3]]=0.707107 (0.707107- 0.707107 b) + 0.707107 (-0.707107 + 
> 0.707107 b)
> 
> Chop works almost as
> Chop[T.M.T]=SparseArray[<5>, {3, 3}]
> But
> Chop[T.M.T][[3,2]]=0
> However
> Chop[T.M.T][[2,3]]=0.707107 (0.707107 - 0.707107 b) + 0.707107 (-0.707107 \
+ 
> 0.707107 b)
> 
> Simplify does nothing. i.e.
> Simplify[T.M.T]=T.M.T
> 
> Chop[Simplify[Normal[T.M.T]]]
> does what I want, but does not work with sparse arrays and in between all \

> elements become large expressions so I run out of memory.
> 
> Maurits Haverkort
> 
> (P.S. Mathematica 5.2 under Linux (Redhat) 1 Gb and Windows XP, 512Mb, the \

> matrix size I aim for is 646*646, however I curently can't evaluate a 
> 134*134 matrix)
> 
> ----- Original Message ----- 
===
> Subject: Re: Zero times a variable in numerical calculations
> 
> 
> 
>> 2) What is the fastest way to multiply a numerical matrix times a
>> half numerical half analytical sparse matrix, whereby the result is
>> band diagonal. (The current result has 0. a + 0. b + 0. c + 0. d +
>> 0. e + 0. f on all places where it should be zero.)
> 
> How are you defining the sparse matrix?
> 
> If I do
> 
> In[17]:=
> a=SparseArray[{{1, 2} -> Random[], {3, 4} -> Random[]},
>   {5, 5}];
> 
> followed by
> 
> In[18]:=
> Normal@a
> 
> Out[18]=
> {{0, 0.49847100218212825, 0, 0, 0}, {0, 0, 0, 0, 0},
>   {0, 0, 0, 0.009744974537243703, 0}, {0, 0, 0, 0, 0},
>   {0, 0, 0, 0, 0}}
> 
> 
> I see the zero elements are exact 0's not inexact 0's. Which means the 
> matrix multiplication will be evaluating 0 a etc which will evaluate to \
0.
> --
> To reply via email subtract one hundred and four 
> 
Hi Maurits,
you might want to manipulate the array rules:
m = SparseArray[{{1, 1} -> a, {2, 2} -> 1., {3, 3} -> 1., {2, 3} -> b, {3, \
2} -> b, {_, _} -> 0}, {3, 3}];
t = SparseArray[{{1, 1} -> 1., {2, 2} -> -Sqrt[0.5], {3, 3} -> Sqrt[0.5], \
{2, 3} -> Sqrt[0.5], {3, 2} -> Sqrt[0.5], {_, _} -> 0}, {3, 3}];
(t . m . t)[[2,3]]
gives the unsimplified element of the matrix.
result = SparseArray[ArrayRules[t . m . t] /.
  HoldPattern[ix_ -> val_] :> ix -> Chop[Simplify[val]]];
result[[2,3]]
gives 0
Normal[result]
-->
{{1.*a, 0, 0}, {0, 1.0000000000000002 - 1.0000000000000002*b, 0},
   {0, 0, 1.0000000000000002 + 1.0000000000000002*b}}
hth,
   Peter
===
Subject: AuthorTools/MakeContents
Hallo,
I have problems with creating a TOC.
Executing MakeContents[nb,\Simple\] aborts with the warning
NotebookSaveWarning::nbmod: The notebook must be saved to use the function 
\\
MakeContents.
although I have saved my notebook before executing the command.
Is it possible to insert a TOC into the notebook that executes this 
command, i.e to insert a notebook into itself.
If not, how is this thing done in every days work. It would be strange if 
notebook A must create a TOC.nb and then try to insert TOC.nb into A.
Is this really the way things are done with MakeContents?
-- 
Oliver Friedrich
My email has no x!
===
Subject: Re: Precision of arguments to FunctionInterpolation
> Here is some further information that makes the strangeness more
> obvious:
> 
> The following three different calls to FunctionInterpolation all fail:
> 
>  FunctionInterpolation[lowPrecisionSin[x], {x, 0, 1.`14}]
>  FunctionInterpolation[lowPrecisionSin[x], {x, 0, 1.`20}]
>  FunctionInterpolation[lowPrecisionSin[x], {x, 0, 1}]
the latter works in Mathematica 5.1 for Windows
> 
> The only way that I can successfully interpolate to an upper bound of
> roughly 1 is by entering the upper bound as the
> exactly-machine-precision number \1.\:
> 
>   FunctionInterpolation[lowPrecisionSin[x], {x, 0, 1.}]
> 
> Can anyone explain why higher-than-machine-precision (e.g. 1.`20),
> lower-than-machine-precision (e.g. 1.`14), and infinite precision (1)
> all fail, but exactly machine-precision succeeds here?
>
No, sorry.
> Secondly, does anyone know the meaning of the InterpolationPrecision
> option to FunctionInterpolation? In particular, with the default
> setting of InterpolationPrecision->Automatic, what value is
> automatically chosen?
> 
> 
You can test this using
Reap[FunctionInterpolation[((Sow[Precision[#1]]; #1) & \
)[lowPrecisionSin[x]], {x, 0., 1}]][[2, 1]]
which gives:
{5.1924023244417254, Infinity, Infinity, 4.999999999999999, \
5.000000000000001, 5., 5., 5.,
   5.000000000000001, 4.999999999999999, 5., 5.000000000000001, 5., 5., \
4.999999999999999, 5.,
   5., 4.999999999999999, 5.000000000000001, 5.000000000000001, \
5.000000000000001, 5., 5.}
If you know the maximum precision of the function you want to interpolate (5 \
in this case), you can use
InterpolationPrecision->5:
intfuns = {f14, f20, finf} = (FunctionInterpolation[lowPrecisionSin[x],
        {x, 0, SetPrecision[1., #1]}, InterpolationPrecision -> 5] & ) /@ \
{14, 20, Infinity};
SameQ @@ intfuns
--> True
Through[{f14, lowPrecisionSin, Sin}[0.6`14.000000000000002]]
--> {0.5646420701356566528`4.127871245474387, \
0.5646424733950353572`5.000000000000002, \
0.5646424733950353572009454457`14.056991706843485}
   Peter
===
Subject: Wolfram Workbench and package development
I've installed and taken a look at Wolfram Workbench (WW).  Before deciding \

to adopt it for my work,  I'd like to raise a few questions concerning how 
one might use it for package development.   It would be interesting to hear \

other people's thoughts.
1. WW is clearly an alternative to the traditional workflow for package 
development in a notebook (which I learned from Roman Maeders' book 
\Programming in Mathematica\). In that scheme, package functions are \
created 
in initialisation cells and the notebook can hold documentation, examples 
and so on.  I've been doing this for years but somehow never got round to 
the next step which is to write additional documentation notebooks for 
integration in the Help Browser.   Presumably this can be done in WW, 
perhaps more efficiently.   It may seem fairly obvious but it would be 
useful if someone would spell out an authoritative recipe for doing so (and \

ideally put it into the help for WW).   Which files should you finally 
\export\ from the project and give to people to install in their 
$(User)BaseDirectory Applications folder?
2. Is the above still the right way to put together and provide Help Browser \

documentation for application packages, even for future versions of 
Mathematica ?
3. What about packages that import other packages, contexts and 
sub-contexts, etc.?
4. When I try to use the Mathematica Import Wizard to bring existing 
notebook-generated packages into WW, I tend to end up with a lot of 
non-InputForm syntax (\\\[Rule], \\\!\\(\, \\\<\,\\\[InvisibleSpace]\ \
etc.).   It 
would be nice to automatically fix/delete these (although perhaps it should \

be the Front End's job).
5. Is it a good or a bad idea to break a single large package over several 
.m files in WW?  If good, is it just a mattering of putting Get[...] 
commands between the BeginPackage and EndPackage ?
6. In debugging or improving a function from an existing, loaded package, I \

often find myself working on a slightly-renamed version that has been copied \

to the Global context (this avoids repeatedly re-launching the kernel, 
particularly when it takes some time to get to the point where the new 
version of the function can be tested).  This method is a little messy, 
especially when it uses private functions.   How would I do this in WW ?
John Jowett
===
Subject: Re:  Wolfram Workbench and package development
I saw your message on the group this morning. I waited for it to come
in via my email account. This evening, I found it in my spam folder...
1)
People have probably mentioned this to you by now, but just in case:
Help files are (supposedly) authored using the author tools palette
and the Help stylesheet. You add help to the browser by using front
end commands in a special file:
one of:
$UserBaseDirectory, $BaseDirectory, $InstallationDirectory
appended with:
\\Applications\\RandomAppName\\Documentation\\English(or your
language)\\BrowserCategories.m
Where BrowserCategories.m contains nested BrowserCategory front end
functions. The syntax of this command may be looked up in the help
browser. If you search the MathGroup for this function name or the
BrowserCategories.m file name, you will find more information.
As noted above, Help files are actually notebooks. Therefore, you
can't edit them directly in Workbench. However, you can certainly
manage the help files as part of a Project in Workbench.
The part of the project that you give to users is the subfolder with
the application name - tell them to put it in (one of:
$UserBaseDirectory, $BaseDirectory,
$InstallationDirectory)\\Applications\\.
You may be wondering what will happen to the init.m file in
YourWorkspace\\YourProject\\AppName\\Kernel if you copy the AppName
folder to your Mathematica Applications directory. The answer is as
follows: Nothing will happen unless you execute <<AppName`. When you
do that, the init.m file will be executed (it usually just calls the
main AppName.m file).
2)
Who knows what they'll do for future version of Mathematica? :-]
3)
Packages may import other packages in the standard way, by appending
the other package contexts to the BeginPackage function call.
4)
Yes, I have noticed this problem as well. A workaround is to use the
Front End's initialization cells autosave package. Mark all the cells
you want as initialization cells, then save. You will get a version of
the text that is pretty good. Alternatively you could cut and paste
using the \plain text\ option. Neither of these options are ideal.
Unless someone knows a better way, we may have to wait for WRI to
produce a better importer.
5)
As you say, it is exactly a matter of putting Get commands in the
appropriate place. No context switching takes place unless you issue
the appropriate commands, AFAIK. I have tested this.
6)
In Workbench, doing that would be prohibitively difficult and also
unnecessary if you learn to use the Debug perspective. In fact, it's
unnecessary in Mathematica as well, if you use the Dialog (and
related) command(s).
The debug perspective as seen when launching a notebook in Debug mode
from Workbench:
The debug subwindow allows you to see the stack (similar to Stack[]
executed in a Mathematica Dialog). The variables subwindow allows you
to see (context free!) symbols that are inside your scope (ie,
variables in a Module statement). There is also an arbitrary
Expressions subwindow ( the help on this window may be out of date).
You can set breakpoints on messages and on arbitrary sections of code
so that you get a chance to actually use those aforementioned
subwindows. You can \step into\ arbitrary expressions (think of adding
to the Stack or going to a sublist in Trace's output). You can also
\step over\ expressions (think of evaluating the next part of a larger
expression like CompoundExpression[a=2+z,Solve[a==x-r,r]] where
Solve[a==x-r,r] is the next part [literally] or moving to the next
item in a list of Trace's output).
I've seen you talk about some pretty advanced stuff on the list, so I
am sure I oversimplified/overexplained.
Anyway, I hope I helped,
> I've installed and taken a look at Wolfram Workbench (WW).  Before \
deciding
> to adopt it for my work,  I'd like to raise a few questions concerning \
how
> one might use it for package development.   It would be interesting to \
hear
> other people's thoughts.
 1. WW is clearly an alternative to the traditional workflow for package
> development in a notebook (which I learned from Roman Maeders' book
> \Programming in Mathematica\). In that scheme, package functions are \
created
> in initialisation cells and the notebook can hold documentation, examples
> and so on.  I've been doing this for years but somehow never got round to
> the next step which is to write additional documentation notebooks for
> integration in the Help Browser.   Presumably this can be done in WW,
> perhaps more efficiently.   It may seem fairly obvious but it would be
> useful if someone would spell out an authoritative recipe for doing so \
(and
> ideally put it into the help for WW).   Which files should you finally
> \export\ from the project and give to people to install in their
> $(User)BaseDirectory Applications folder?
 2. Is the above still the right way to put together and provide Help \
Browser
> documentation for application packages, even for future versions of
> Mathematica ?
 3. What about packages that import other packages, contexts and
> sub-contexts, etc.?
 4. When I try to use the Mathematica Import Wizard to bring existing
> notebook-generated packages into WW, I tend to end up with a lot of
> non-InputForm syntax (\\\[Rule], \\\!\\(\, \
\\\<\,\\\[InvisibleSpace]\ etc.).   It
> would be nice to automatically fix/delete these (although perhaps it \
should
> be the Front End's job).
 5. Is it a good or a bad idea to break a single large package over \
several
> .m files in WW?  If good, is it just a mattering of putting Get[...]
> commands between the BeginPackage and EndPackage ?
 6. In debugging or improving a function from an existing, loaded package, \
I
> often find myself working on a slightly-renamed version that has been \
copied
> to the Global context (this avoids repeatedly re-launching the kernel,
> particularly when it takes some time to get to the point where the new
> version of the function can be tested).  This method is a little messy,
> especially when it uses private functions.   How would I do this in WW ?
 John Jowett


>
-- 
http://chris.chiasson.name/
===
Subject: Re: Wolfram Workbench and package development
> I've installed and taken a look at Wolfram Workbench (WW).  Before \
deciding
> to adopt it for my work,  I'd like to raise a few questions concerning \
how
> one might use it for package development.   It would be interesting to \
hear
> other people's thoughts.
 1. WW is clearly an alternative to the traditional workflow for package
> development in a notebook (which I learned from Roman Maeders' book
> \Programming in Mathematica\). In that scheme, package functions are \
created
> in initialisation cells and the notebook can hold documentation, examples
> and so on.  I've been doing this for years but somehow never got round to
> the next step which is to write additional documentation notebooks for
> integration in the Help Browser.   Presumably this can be done in WW,
> perhaps more efficiently.   It may seem fairly obvious but it would be
> useful if someone would spell out an authoritative recipe for doing so \
(and
> ideally put it into the help for WW).   Which files should you finally
> \export\ from the project and give to people to install in their
> $(User)BaseDirectory Applications folder?
I am not sure if you have looked at the help browser, but these topics
are well covered there.  integration in the Help browser is more of a
Mathematica related issue.
 > 2. Is the above still the right way to put together and provide Help
Browser
> documentation for application packages, even for future versions of
> Mathematica ?
 3. What about packages that import other packages, contexts and
> sub-contexts, etc.?
 4. When I try to use the Mathematica Import Wizard to bring existing
> notebook-generated packages into WW, I tend to end up with a lot of
> non-InputForm syntax (\\\[Rule], \\\!\\(\, \
\\\<\,\\\[InvisibleSpace]\ etc.).   It
> would be nice to automatically fix/delete these (although perhaps it \
should
> be the Front End's job).
> 5. Is it a good or a bad idea to break a single large package over \
several
> .m files in WW?  If good, is it just a mattering of putting Get[...]
> commands between the BeginPackage and EndPackage ?
I am not a developer for Workbench, the idea here is that Workbench
essentially creates  .m files, so all the techniques and practises that
one uses to create packages in Mathematica apply here as well.
> 6. In debugging or improving a function from an existing, loaded package, \
I
> often find myself working on a slightly-renamed version that has been \
copied
> to the Global context (this avoids repeatedly re-launching the kernel,
> particularly when it takes some time to get to the point where the new
> version of the function can be tested).  This method is a little messy,
> especially when it uses private functions.   How would I do this in WW ?
If you just save your work in Workbench, this gets automatically
updated to your Mathematica session. If you then Run/Debug your
file--these changes are automatically reflected
Hope this helps
Pratik 
Wolfram Research 
> 
> John Jowett
===
Subject: Change CellTags with FrontEnd
I'd like to change the celltags of several cells that are not adjacent
via a program (not via the menu). Doing this for one cell is no
problem: I use a rule like
{Cell[a_, st_String, b___, CellTags -> tgs_, c___] :>
    Cell[a, st, b, CellTags -> Union[Flatten[{tgs, newtag}]], c],
Cell[a_, st_String, b___ /; FreeQ[{b}, CellTags, \\[Infinity]]] :>
    Cell[a, st, b, CellTags -> settag]}.
The problem arises now if several cells are selected: I read them in,
apply the rule and paste them back, but the new cells appear now at the
wrong place. One can simply imitate this behaviour of M by selecting
several not adjacent cells and doing copy and paste: One will see that
the cells are pasted at the wrong place.
Doing the same thing via the menu entry Find > Add/Remove Cell Tags
this is no problem at all. So it appears to me that a solution using
the FrontEnd would be the best thing. But: I can find no documentation
about a FrontEndToken like \AddCellTags\. The best thing I can do is
use
FrontEndTokenExecute[CellTagsEditDialog]
to open the dialog to change CellTags, but this is not what I want.
So can anybody tell me if there is an additional FrontEndToken that can
do what I want or has another idea? It should work for M4.0 - M5.2.
===
Subject: Re: Zero times a variable in numerical calculations
>Chop[0.0 a] indeed works and returns nicely 0 however I still have
>some problems when using this on Sparse Arrays.
>I define: M = SparseArray[{{1, 1} -> a, {2, 2} -> 1.0, {3, 3} -1.0, {2, 3} -> b, {3, 2} -> b, {_, _} -> 0}, {3, 3}];
You don't need the {_,_}->0 part since by default SparseArray assumes all \
elements not given specific values are 0. However, deleting this won't solve \
your problem.
<snipChop[Simplify[Normal[T.M.T]]] does what I want, but does not work
>with sparse arrays and in between all elements become large
>expressions so I run out of memory.
You don't have to convert the SparseArray to normal form in order to use \
Chop and Simplify. You could use pattern matching to operate directly on the \
defined elements, i.e.,
In[3]:=
mat = T . M . T /. SparseArray[a_, b_, c_, {x__, y_}] :> 
     SparseArray[a, b, c, {x, Chop[Simplify[y]]}]; 
Normal[mat]
Out[4]=
{{1.*a, 0, 0}, {0, 1.0000000000000002 - 1.0000000000000002*b, 0}, 
  {0, 0, 1.0000000000000002*b + 1.0000000000000002}}
--
To reply via email subtract one hundred and four
===
Subject: Re:  Programming with .NET / hand over of variables
i don't know how simple or difficult this question is, but i do not know
>how to work this problem out.
>i have an problem with the implementation of Mathematica in my
>.NET-Program.
Here is the code:
mathKernel.Compute(\puffer1 =
>ReadList[\\\v3504b.txt\\\,{Number,Number}]\);  //Read an array of \
numbers
>from a file
>mathErgebnis = new
>Wolfram.NETLink.Expr(mathKernel.Result);                       //The
>Result, which is send back by mathKernel is saved as an Expression
>double[][] ar = ?
>                                           //Transform the Result in an
>array which i can handle with C#
I want to have the 2*N array, which is saved in puffer1, transformed in
>an C# equivalent.
>What have i to do, so that it works?
Andreas,
Since you want the result as an Expr, set the ResultFormat property to 
ensure that:
     mathKernel.ResultFormat = MathKernel.ResultFormatType.Expr;
     mathKernel.Compute(\puffer1 = \
ReadList[\\\v3504b.txt\\\,{Number,Number}]\);
     mathErgebnis = (Expr) mathKernel.Result;
To get the Expr's value as a native C# array, use the AsArray() method:
     double[,] ar = mathErgebnis.AsArray(ExpressionType.Real, 2);
You say you want a double[][] instead of a double[,]. If you really want 
the \array of arrays\ format you will have to manually allocate a 
double[][] of the correct size and fill it from the elements in ar. It is 
generally more convenient to work with multidimensional arrays (double[,]) 
than arrays-of-arrays (double[][]) in C#.
Todd Gayley
Wolfram Research
===
Subject: Re:  mathematica 4.0 + NET/Link doesn't work?
actually i have mathematica in version 4.0 installed. Now i want to
>build a UI for a given mathematica program.
>I have downloaded MathLink, J/Link and NET/Link from www.wolfram.com,
>and so they should run in the latest version.
>But when i want to initialize the Link, i got an error:
>Here is the code:
>---------------------------------------------------------------------------\
---------------------------------
>       In[1]:= <<NETLink`
       In[2]:= InstallNET[]
LinkOpen::\linke\: \MathLink error: \\!\\(\\\the linkname was missing \
or not in \\
>the proper form\\\\\)\
InstallNET::\fail\: \A link to the .NET runtime could not be \
established.
>---------------------------------------------------------------------------\
-----------------------------------
>i become desperate, i need that program to work.
Andreas,
Here's a quick one-line fix to get .NET/Link to work back in Mathematica \
4.0:
     SetOptions[LinkOpen, LinkMode->Launch]
Now InstallNET[] will succeed. I can't guarantee that everything in 
.NET/Link will work in Mathematica 4.0, but then neither can I think of 
anything that will fail.
Todd Gayley
Wolfram Research 
===
Subject: Something in the spirit of letrec
Hi all,
I'm trying to somehow come up with a function WithSelf such that
evaluating
  WithSelf[{a,self[1]+b,self[2]+2}]
will result in {a,a+b,a+b+c}. The idea is that the each
item in the argument may depend on the previous ones;
so for example \self[1]+b\ means \what's in the first position
in this list, plus b\. There will be no \forward references\
in the argument.
I've tried with
  SetAttributes[WithSelf, HoldFirst];
  WithSelf[l_, self_] := Module[{ll, r = {}},
      ll = Function[self, #] & /@ l;
      Scan[AppendTo[r, #[r]] &, ll];
      r
      ]
This gives the correct correct results, but emits warnings:
for example,
  WithSelf[{1, self[[1]] + 2, self[[2]] + 3, self[[1]]}, self]
results in {1, 3, 6, 1}, which is OK, but emits a few of
  Part::partd: Part specification self[[1]] is longer than depth of
object.
I guess there is some evaluation going on which I am
not seeing...
I'm quite sure what I am trying to do is possible... but how?
Ideas?
-- m
===
Subject: How to do Chebyshev expansion in Mathematica?
===
Subject: near Planck's mass
Back about 1987 or before
 I was studying gravitation and quantum mechanics
and I came across  this:
(* Classical electromagnetic radius*)
r = e^2/(m*c^2)
(* Riemannian Curvature*)
R = -2/r^2
(* 3 space volume*)
V = (4*Pi/3)*r^3
(*scalar energy densidty*)
T = m*c^2/V
(* solution for R = -8Pi*G*T(v, v) : scalar reduction of Einstein's
        general relativity*)
Solve[R + (8*Pi*G/c^4)*T == 0, m]
e = 4.80325*10^(-10)
G = 6.6732*10^(-8)
(* mass in grams*)
m = e/Sqrt[3*G]
1.073513458510097`*10^-6
At the time I hadn't run into Planck's mass.
If you put in
 r=h/(m*c)
in this calculation  you get a version of the Planck mass.
I was wondering if this was a well known calculation like Planck's mass?
Roger
===
Subject: Re:  returning a variable's name, rather than the variable's \
contents
 There must be a simple way to do this but it eludes me. Take for  
> example the
> following:

 In[2]:= a=1;b=2;c=3;

 In[3]:= Max[a,b,c]

 Out[3]= 3

 What would I do if I wanted Out[3] to equal \c\ ?



> Michael Stern

Here is a real newbie approach without assignments:
In[1]:=
eqns = {a -> 1, b -> 2,
    c -> 3}
revs = {1 -> a, 2 -> b,
    3 -> c}
Out[1]=
{a -> 1, b -> 2, c -> 3}
Out[2]=
{1 -> a, 2 -> b, 3 -> c}
In[3]:=
Max[a, b, c] /. eqns /. revs
Out[3]=
c
J\.87nos
===
Subject: Re:  returning a variable's name, rather than the variable's \
contents
variables = {a, b, c, d, e};
values = {1, 2, 3, 4, 5};
Max[values] /. Thread[values -> variables]
> There must be a simple way to do this but it eludes me. Take for example \
the
> following:
> 
>  
> 
> In[2]:= a=1;b=2;c=3;
> 
>  
> 
> In[3]:= Max[a,b,c]
> 
>  
> 
> Out[3]= 3
> 
>  
> 
> What would I do if I wanted Out[3] to equal \c\ ?
> 
>  
> 
> 
>  
> 
> Michael Stern
> 
> 
> 
-- 
Takashi Yoshino
===
Subject: returning a variable's name, rather than the variable's contents
There must be a simple way to do this but it eludes me. Take for example \
the
following:
 
In[2]:= a=1;b=2;c=3;
 
In[3]:= Max[a,b,c]
 
Out[3]= 3
 
What would I do if I wanted Out[3] to equal \c\ ?
 
 
Michael Stern
===
Subject: Re: returning a variable's name, rather than the variable's \
contents
> There must be a simple way to do this but it eludes me. Take for example \
the
> following: 
> 
> In[2]:= a=1;b=2;c=3;
> 
> In[3]:= Max[a,b,c]
> 
> Out[3]= 3
> 
> What would I do if I wanted Out[3] to equal \c\ ?
What about the following trick?
In[1]:=
vars = {a, b, c};
varNames = ToString /@ vars;
values = {1, 2, 3};
vars = values;
First[Extract[varNames, Position[vars, Max[vars]]]]
Out[5]=
\c\
HTH,
Jean-Marc
===
Subject: Re:  Converting Mathematics slides into PDF
 I read some of the posts on this topic, but I don't think they  
> address the
> problem I am having.
 Is there a way to generate PDF file of Mathematica slides such that  
> what
> shows on the screen occupies one page or less in the PDF file. For  
> example,
> I have graphs which do not even fill the screen completely, but do  
> not fit
> in the page when I make the PDF file (that is, only a part shows up  
> on the
> page and other part is cut out). Similarly, some of the text that  
> fits in
> half the screen fills up multiple pages in the PDF file.

> Bharat.
I think when you create a pdf it uses the current printer driver for  
it.  So make sure that you have the right printer in mind and page  
setup is correct.
J\.87nos
===
Subject: Converting Mathematics slides into PDF
I read some of the posts on this topic, but I don't think they address the
problem I am having.
Is there a way to generate PDF file of Mathematica slides such that what
shows on the screen occupies one page or less in the PDF file. For example,
I have graphs which do not even fill the screen completely, but do not fit
in the page when I make the PDF file (that is, only a part shows up on the
page and other part is cut out). Similarly, some of the text that fits in
half the screen fills up multiple pages in the PDF file.
Bharat.
===
Subject: Re: Converting Mathematics slides into PDF
> I read some of the posts on this topic, but I don't think they address \
the
> problem I am having.
 Is there a way to generate PDF file of Mathematica slides such that what
> shows on the screen occupies one page or less in the PDF file. For \
example,
> I have graphs which do not even fill the screen completely, but do not \
fit
> in the page when I make the PDF file (that is, only a part shows up on \
the
> page and other part is cut out). Similarly, some of the text that fits in
> half the screen fills up multiple pages in the PDF file.
What platform are you using?  How are you generating the PDF file?  Are
you using Mathematica's built in Export[] command?  Mac OS X's Print
Preview?  Acrobat print driver?
-Rob
===
Subject: Re:  Re: Pattern match question
----- Original Message ----- 
===
Subject:  Re: Pattern match question
> hi
>> txt = \ZACCZBNRCSAACXBXX\;
>> letters = \ABC\;
>> i want to find the first occurrences of any of the
>> six combinations of the letters of the set \ABC\ Globally, and
>> without overlap option. and the space between letters does not
>> important.
>> in the above txt string the result must be:
>> Out[]:=
>> ACCZB
>> CSAACXB
>> i wish a solution using mathematica regular expressions.
>> the Regex pattern  (A|B|C).*?(A|B|C).*?(A|B|C)  will give the out:
>> ACC ,  BNRCSA  ,  ACXB   because it considers the permutations
>> and not the combinations
>> the following is an old fashion program which will emulate the human
>> pencil and
>> paper method,  will solve the problem, but i am sure there are a better
>> solutions.
>> txt = \ZACCZBNRCSAACXBXX\;
>> letters = \ABC\;
>> ptrnLtrs = \\;
>> (* make a string of 26 zero's as the number of the alphbet*)
>> For[i = 1, i <= 26, ptrnLtrs = StringJoin[ptrnLtrs, \0\]; i++]
>> (* replace every letter of the pattern letters *)
>> (* with a corresponding 1 in the string of the zero's *)
>> For[i = 1, i <= StringLength[letters],
>>     num = ToCharacterCode[StringTake[letters, {i, i}]];
>>     num = num - 64;
>>     ptrnLtrs = StringReplacePart[ptrnLtrs, \1\, Flatten[{num, num}]];
>>     i++];
>> (* the procedural pattern match *)
>> ptrnLtrsBak = ptrnLtrs; y = 0;    (* backup for the ptrnLtrs *)
>> beginFlag = 0; result = \\; lst = {};
>> For[i = 1, i <= StringLength[txt],
>>     OneLetter = StringTake[txt, {i, i}];
>>     If[beginFlag == 0 && StringCases[letters, OneLetter] == {},
>>        Goto[jmp]];
>>      num = ToCharacterCode[StringTake[txt, {i, i}]] - 64;
>>      If[StringTake[ptrnLtrs, num] == \1\,
>>         result = StringJoin[result, OneLetter];
>>         ptrnLtrs = StringReplacePart[ptrnLtrs, \0\, Flatten[{num,
>> num}]];
>>         , result = StringJoin[result, OneLetter];];
>>       beginFlag = 1;
>>       If[ToExpression[ptrnLtrs] == 0, ptrnLtrs = ptrnLtrsBak;
>>           Print[result];
>>           result = \\;  beginFlag = 0;];
>>        Label[jmp];
>>        i++]
>> Out[]:=
>> ACCZB
>> CSAACXB
>> peter glassy
 Here's a solution that uses string expressions:
 In[1]:= Module[
>  {Lpatt = StringExpression @@@ (Insert[#, ___, {{2}, {3}}]&) /@
>      Permutations@ {\A\, \B\, \C\}},
>  StringCases[\ZACCZBNRCSAACXBXX\,
>    ShortestMatch[s__] /; StringMatchQ[s, Lpatt]]]
 Out[1]= {\ACCZB\, \CSAACXB\}
 Note that in StringCases a list of patterns {patt1, patt2, ...} is
> equivalent to patt1 | patt2 | ... . We cannot directly use
> ShortestMatch[patt1 | patt2] because this merely makes all the
> quantifiers in the regex lazy but doesn't guarantee that we get the
> shortest possible match:
 In[2]:= StringPattern`PatternConvert[
>  ShortestMatch[(\A\ ~~ ___ ~~ \B\) | (\A\ ~~ ___ ~~ \C\)]]
 Out[2]= {\(?ms)(?:A.*?B|A.*?C)\, {}, {}, Hold[None]}
 In[3]:= StringCases[\ACB\,
>  ShortestMatch[(\A\ ~~ ___ ~~ \B\) | (\A\ ~~ ___ ~~ \C\)]]
 Out[3]= {\ACB\}
 The shortest match would be \AC\. So it's interesting to consider how
> we can obtain the same answer {\ACCZB\, \CSAACXB\} using only
> RegularExpression without external conditions.
 Maxim Rytin
> m.r@inbox.ru
Yes, because of the underlying algorithm used, it's true that ShortestMatch \

does not affect Alternatives, only quantifiers like \__\ and \..\.
Here are a couple of regular expression approaches:
In[1]:= StringCases[\ZACCZBNRCSAACXBXX\,
RegularExpression[
\A.*?(?:B.*?C|C.*?B)|B.*?(?:A.*?C|C.*?A)|C.*?(?:A.*?B|B.*?A)\]]
Out[1]= {\ACCZB\, \CSAACXB\}
In[2]:= StringCases[\ZACCZBNRCSAACXBXX\,
RegularExpression[\(A|B|C).*?(?!\\\\1)(A|B|C).*?(?!\\\\1)(?!\\\\2)(?:A|B|C)\
\]]
Out[2]= {\ACCZB\, \CSAACXB\}
The first is similar to Maxim's pattern, except it makes sure that each 
alternative can only be matched once, this also helps with efficiency. It 
can also be written as a StringExpression (see below).
The second is similar to Peter's original pattern, but it uses a \negative \

lookahead\ to make sure that the same character is not repeated. There is \
no 
StringExpression equivalent of this.
The first pattern is probably slightly faster in general (if it's confusing \

to read, remember that .*? is basically ShortestMatch[___] and (?:...) is 
the way to group things without wasting time on storing the subpattern). 
Here is the StringExpression version of it:
In[3]:=
StringCases[\ZACCZBNRCSAACXBXX\,
ShortestMatch[
(\A\ ~~ ___ ~~ ((\B\ ~~ ___ ~~ \C\) | (\C\ ~~ ___ ~~ \B\))) |
(\B\ ~~ ___ ~~ ((\A\ ~~ ___ ~~ \C\) | (\C\ ~~ ___ ~~ \A\))) |
(\C\ ~~ ___ ~~ ((\A\ ~~ ___ ~~ \B\) | (\B\ ~~ ___ ~~ \A\)))]]
Out[3]= {\ACCZB\, \CSAACXB\}
Oyvind Tafjord
Wolfram Research 
===
Subject: Re: Pattern match question
> hi
> txt = \ZACCZBNRCSAACXBXX\;
> letters = \ABC\;
> i want to find the first occurrences of any of the
> six combinations of the letters of the set \ABC\ Globally, and
> without overlap option. and the space between letters does not
> important.
> in the above txt string the result must be:
> Out[]:=
> ACCZB
> CSAACXB
> i wish a solution using mathematica regular expressions.
> the Regex pattern  (A|B|C).*?(A|B|C).*?(A|B|C)  will give the out:
> ACC ,  BNRCSA  ,  ACXB   because it considers the permutations
> and not the combinations
> the following is an old fashion program which will emulate the human
> pencil and
> paper method,  will solve the problem, but i am sure there are a better
> solutions.
 txt = \ZACCZBNRCSAACXBXX\;
> letters = \ABC\;
> ptrnLtrs = \\;
> (* make a string of 26 zero's as the number of the alphbet*)
> For[i = 1, i <= 26, ptrnLtrs = StringJoin[ptrnLtrs, \0\]; i++]
> (* replace every letter of the pattern letters *)
> (* with a corresponding 1 in the string of the zero's *)
> For[i = 1, i <= StringLength[letters],
>     num = ToCharacterCode[StringTake[letters, {i, i}]];
>     num = num - 64;
>     ptrnLtrs = StringReplacePart[ptrnLtrs, \1\, Flatten[{num, num}]];
>     i++];
 (* the procedural pattern match *)
> ptrnLtrsBak = ptrnLtrs; y = 0;    (* backup for the ptrnLtrs *)
> beginFlag = 0; result = \\; lst = {};
> For[i = 1, i <= StringLength[txt],
>     OneLetter = StringTake[txt, {i, i}];
>     If[beginFlag == 0 && StringCases[letters, OneLetter] == {},
>        Goto[jmp]];
>      num = ToCharacterCode[StringTake[txt, {i, i}]] - 64;
>      If[StringTake[ptrnLtrs, num] == \1\,
>         result = StringJoin[result, OneLetter];
>         ptrnLtrs = StringReplacePart[ptrnLtrs, \0\, Flatten[{num,
> num}]];
>         , result = StringJoin[result, OneLetter];];
>       beginFlag = 1;
>       If[ToExpression[ptrnLtrs] == 0, ptrnLtrs = ptrnLtrsBak;
>           Print[result];
>           result = \\;  beginFlag = 0;];
>        Label[jmp];
>        i++]
 Out[]:=
> ACCZB
> CSAACXB
 peter glassy
Here's a solution that uses string expressions:
In[1]:= Module[
  {Lpatt = StringExpression @@@ (Insert[#, ___, {{2}, {3}}]&) /@
      Permutations@ {\A\, \B\, \C\}},
  StringCases[\ZACCZBNRCSAACXBXX\,
    ShortestMatch[s__] /; StringMatchQ[s, Lpatt]]]
Out[1]= {\ACCZB\, \CSAACXB\}
Note that in StringCases a list of patterns {patt1, patt2, ...} is
equivalent to patt1 | patt2 | ... . We cannot directly use
ShortestMatch[patt1 | patt2] because this merely makes all the
quantifiers in the regex lazy but doesn't guarantee that we get the
shortest possible match:
In[2]:= StringPattern`PatternConvert[
  ShortestMatch[(\A\ ~~ ___ ~~ \B\) | (\A\ ~~ ___ ~~ \C\)]]
Out[2]= {\(?ms)(?:A.*?B|A.*?C)\, {}, {}, Hold[None]}
In[3]:= StringCases[\ACB\,
  ShortestMatch[(\A\ ~~ ___ ~~ \B\) | (\A\ ~~ ___ ~~ \C\)]]
Out[3]= {\ACB\}
The shortest match would be \AC\. So it's interesting to consider how
we can obtain the same answer {\ACCZB\, \CSAACXB\} using only
RegularExpression without external conditions.
Maxim Rytin
m.r@inbox.ru
===
Subject: Re: Pattern match question
> hi
> txt = \ZACCZBNRCSAACXBXX\;
> letters = \ABC\;
> i want to find the first occurrences of any of the
> six combinations of the letters of the set \ABC\ Globally, and
> without overlap option. and the space between letters does not
> important.
> in the above txt string the result must be:
> Out[]:=
> ACCZB
> CSAACXB
> i wish a solution using mathematica regular expressions.
> the Regex pattern  (A|B|C).*?(A|B|C).*?(A|B|C)  will give the out:
> ACC ,  BNRCSA  ,  ACXB   because it considers the permutations
> and not the combinations
> the following is an old fashion program which will emulate the human
> pencil and
> paper method,  will solve the problem, but i am sure there are a better
> solutions.
 txt = \ZACCZBNRCSAACXBXX\;
> letters = \ABC\;
> ptrnLtrs = \\;
> (* make a string of 26 zero's as the number of the alphbet*)
> For[i = 1, i <= 26, ptrnLtrs = StringJoin[ptrnLtrs, \0\]; i++]
> (* replace every letter of the pattern letters *)
> (* with a corresponding 1 in the string of the zero's *)
> For[i = 1, i <= StringLength[letters],
>     num = ToCharacterCode[StringTake[letters, {i, i}]];
>     num = num - 64;
>     ptrnLtrs = StringReplacePart[ptrnLtrs, \1\, Flatten[{num, num}]];
>     i++];
 (* the procedural pattern match *)
> ptrnLtrsBak = ptrnLtrs; y = 0;    (* backup for the ptrnLtrs *)
> beginFlag = 0; result = \\; lst = {};
> For[i = 1, i <= StringLength[txt],
>     OneLetter = StringTake[txt, {i, i}];
>     If[beginFlag == 0 && StringCases[letters, OneLetter] == {},
>        Goto[jmp]];
>      num = ToCharacterCode[StringTake[txt, {i, i}]] - 64;
>      If[StringTake[ptrnLtrs, num] == \1\,
>         result = StringJoin[result, OneLetter];
>         ptrnLtrs = StringReplacePart[ptrnLtrs, \0\, Flatten[{num,
> num}]];
>         , result = StringJoin[result, OneLetter];];
>       beginFlag = 1;
>       If[ToExpression[ptrnLtrs] == 0, ptrnLtrs = ptrnLtrsBak;
>           Print[result];
>           result = \\;  beginFlag = 0;];
>        Label[jmp];
>        i++]
 Out[]:=
> ACCZB
> CSAACXB
 peter glassy
 Here's a solution that uses string expressions:
 In[1]:= Module[
>   {Lpatt = StringExpression @@@ (Insert[#, ___, {{2}, {3}}]&) /@
>       Permutations@ {\A\, \B\, \C\}},
>   StringCases[\ZACCZBNRCSAACXBXX\,
>     ShortestMatch[s__] /; StringMatchQ[s, Lpatt]]]
 Out[1]= {\ACCZB\, \CSAACXB\}
 Note that in StringCases a list of patterns {patt1, patt2, ...} is
> equivalent to patt1 | patt2 | ... . We cannot directly use
> ShortestMatch[patt1 | patt2] because this merely makes all the
> quantifiers in the regex lazy but doesn't guarantee that we get the
> shortest possible match:
 In[2]:= StringPattern`PatternConvert[
>   ShortestMatch[(\A\ ~~ ___ ~~ \B\) | (\A\ ~~ ___ ~~ \C\)]]
 Out[2]= {\(?ms)(?:A.*?B|A.*?C)\, {}, {}, Hold[None]}
 In[3]:= StringCases[\ACB\,
>   ShortestMatch[(\A\ ~~ ___ ~~ \B\) | (\A\ ~~ ___ ~~ \C\)]]
 Out[3]= {\ACB\}
 The shortest match would be \AC\. So it's interesting to consider how
> we can obtain the same answer {\ACCZB\, \CSAACXB\} using only
> RegularExpression without external conditions.
 Maxim Rytin
> m.r@inbox.ru
Here are two shorter solutions:
In[1]:= StringCases[\ZACCZBNRCSAACXBXX\,
  s : ({\A\, \B\, \C\} ~~ ShortestMatch[__]) /;
Out[1]= {\ACCZB\, \CSAACXB\}
In[2]:= StringCases[\ZACCZBNRCSAACXBXX\,
  \
RegularExpression[\(A|B|C).*?(?!\\\\1)(A|B|C).*?(?!\\\\1|\\\\2)(A|B|C)\]]
Out[2]= {\ACCZB\, \CSAACXB\}
Maxim Rytin
m.r@inbox.ru
===
Subject: Re: How to treat this false singular point?
> I have a function:
> f[x_]:=(x-x1)Log[Abs[x-x1]] + (x-x2)Log[Abs[x-x2]] + ... +
> (x-xn)Log[Abs[x-xn]],
> {x1,x2,...,xn}={100,200,300,...} for instance
> How to get value: f[x] as there are different singular at different x?
> I know at x=xn, f[x]==1, But Mathematica return: \Indeterminate\, What
> should I do?
>  what others do in C++, Fortran ?
The easiest way would be
  Unprotect[Log]; Log[0] = 0; Log[0.] = 0; Protect[Log];
That would do what you want here, but might cause problems elsewhere.
Another way would be to define a new function
  xLog[0] = 0; xLog[0.] = 0; xLog[x_] := x*Log[Abs[x]]
Then
  f[x_] := Tr@Map[xLog,x-{100,200,...}]
===
Subject: Re:  How to treat this false singular point?
Use Limit.  For example, define f[x_]
f[x_] := Limit[(xi - x1)Log[Abs[xi -
            x1]] + (xi - x2)Log[Abs[xi - x2]], xi -> x, Assumptions -> {
            Element[x1, Reals], Element[x2, Reals]}]
and estimate f[x1].  The output is like
((x1 - x2)*Log[(x1 - x2)^2])/2
> I have a function:
> f[x_]:=(x-x1)Log[Abs[x-x1]] + (x-x2)Log[Abs[x-x2]] + ... +
> (x-xn)Log[Abs[x-xn]],
> {x1,x2,...,xn}={100,200,300,...} for instance
> How to get value: f[x] as there are different singular at different x?
> I know at x=xn, f[x]==1, But Mathematica return: \Indeterminate\, What
> should I do?
>  what others do in C++, Fortran ?
> 
-- 
Takashi Yoshino
===
Subject: Re: How to treat this false singular point?
> I have a function:
> f[x_]:=(x-x1)Log[Abs[x-x1]] + (x-x2)Log[Abs[x-x2]] + ... +
> (x-xn)Log[Abs[x-xn]],
> {x1,x2,...,xn}={100,200,300,...} for instance
> How to get value: f[x] as there are different singular at different x?
> I know at x=xn, f[x]==1, But Mathematica return: \Indeterminate\, What
> should I do?
>  what others do in C++, Fortran ?
Use some *If* statements around each term for instance.
In[1]:=
f[x_] := If[x != 100, (x - 100)*Log[Abs[x - 100]], 0] +
          If[x != 200, (x - 200)*Log[Abs[x - 200]], 0] +
          If[x != 300, (x - 300)*Log[Abs[x - 300]], 0]
In[2]:=
f[50]
Out[2]=
-50*Log[50] - 150*Log[150] - 250*Log[250]
In[3]:=
f[100.]
Out[3]=
-1520.18
In[4]:=
f[200]
Out[4]=
0
In[5]:=
Plot[f[x], {x, 0, 500}];
HTH,
Jean-Marc
===
Subject: Re:  How to treat this false singular point?
f[x_, xList_List] := Total[(x-xList)*Log[Abs[x-xList]]];
f[x, {x1, x2, x3}]
(x - x1)*Log[Abs[x - x1]] + (x - x2)*Log[Abs[x - x2]] + (x - x3)*Log[Abs[x - \
x3]]
a singularity will exist for each value of x in xList; therefor, f[xn, \
xList] cannot have the value 1
f[5. - 10^-15, Range[5]]
10.2273
Plot[f[x,Range[5]],{x,-1,5}];
Bob Hanlon
> I have a function:
> f[x_]:=(x-x1)Log[Abs[x-x1]] + (x-x2)Log[Abs[x-x2]] + ... +
> (x-xn)Log[Abs[x-xn]],
> {x1,x2,...,xn}={100,200,300,...} for instance
> How to get value: f[x] as there are different singular at different x?
> I know at x=xn, f[x]==1, But Mathematica return: \Indeterminate\, What
> should I do?
>  what others do in C++, Fortran ?
> 
===
Subject: Re: How to treat this false singular point?
I had made a mistake again, haha, sorry.
In: Limit[(x-10)*Log[Abs[x-10]],x->10]
Out: 0
> f[x_, xList_List] := Total[(x-xList)*Log[Abs[x-xList]]];
 f[x, {x1, x2, x3}]
 (x - x1)*Log[Abs[x - x1]] + (x - x2)*Log[Abs[x - x2]] + (x - x3)*Log[Abs[x \
- x3]]
 a singularity will exist for each value of x in xList; therefor, f[xn, \
xList] cannot have the value 1
 f[5. - 10^-15, Range[5]]
 10.2273
 Plot[f[x,Range[5]],{x,-1,5}];

> Bob Hanlon
 I have a function:
> f[x_]:=(x-x1)Log[Abs[x-x1]] + (x-x2)Log[Abs[x-x2]] + ... +
> (x-xn)Log[Abs[x-xn]],
> {x1,x2,...,xn}={100,200,300,...} for instance
> How to get value: f[x] as there are different singular at different x?
> I know at x=xn, f[x]==1, But Mathematica return: \Indeterminate\, \
What
> should I do?
>  what others do in C++, Fortran ?
>
===
Subject: Re: How to treat this false singular point?
> f[x_, xList_List] := Total[(x-xList)*Log[Abs[x-xList]]];
 f[x, {x1, x2, x3}]
 (x - x1)*Log[Abs[x - x1]] + (x - x2)*Log[Abs[x - x2]] + (x - x3)*Log[Abs[x \
- x3]]
 a singularity will exist for each value of x in xList; therefor, f[xn, \
xList] cannot have the value 1
 f[5. - 10^-15, Range[5]]
 10.2273
 Plot[f[x,Range[5]],{x,-1,5}];

> Bob Hanlon
Sorry, I had made a mistake, I mean the limit of (x-xn)Log[x-xn]==1
when x->xn and x>xn.
The f[x] is a part of another function g[x], Mathematica didn't get
this value when met x=xn which I thought it should.
===
Subject: RE:  Change CellTags with FrontEnd
Ingolf Dahl has an extensive package for manipulating cell tags. It is
available at
And a new version of my package is available at
http://web.telia.com/~u31815170/Mathematica/
and his email address for more information is ingolf.dahl@telia.com.
David Park
djmp@earthlink.net
http://home.earthlink.net/~djmp/
I'd like to change the celltags of several cells that are not adjacent
via a program (not via the menu). Doing this for one cell is no
problem: I use a rule like
{Cell[a_, st_String, b___, CellTags -> tgs_, c___] :>
    Cell[a, st, b, CellTags -> Union[Flatten[{tgs, newtag}]], c],
Cell[a_, st_String, b___ /; FreeQ[{b}, CellTags, \\[Infinity]]] :>
    Cell[a, st, b, CellTags -> settag]}.
The problem arises now if several cells are selected: I read them in,
apply the rule and paste them back, but the new cells appear now at the
wrong place. One can simply imitate this behaviour of M by selecting
several not adjacent cells and doing copy and paste: One will see that
the cells are pasted at the wrong place.
Doing the same thing via the menu entry Find > Add/Remove Cell Tags
this is no problem at all. So it appears to me that a solution using
the FrontEnd would be the best thing. But: I can find no documentation
about a FrontEndToken like \AddCellTags\. The best thing I can do is
use
FrontEndTokenExecute[CellTagsEditDialog]
to open the dialog to change CellTags, but this is not what I want.
So can anybody tell me if there is an additional FrontEndToken that can
do what I want or has another idea? It should work for M4.0 - M5.2.
===
Subject: RE:  returning a variable's name, rather than the variable's \
contents
Michael,
I like this question because it involves the question of how to work
symbolically and numerically at the same time. Doing that is often the best
way to use Mathematica.
One answer is not to set values for your symbolic variables but do use \
rules
instead. Then you can have the best of both worlds. Here is one method for
doing that.
data = {a -> 1, b -> 2, c -> 3};
First@First@Sort[data, Part[#1, 2] >= Part[#2, 2] &]
c
David Park
djmp@earthlink.net
http://home.earthlink.net/~djmp/
There must be a simple way to do this but it eludes me. Take for example \
the
following:
In[2]:= a=1;b=2;c=3;
In[3]:= Max[a,b,c]
Out[3]= 3
What would I do if I wanted Out[3] to equal \c\ ?
Michael Stern
===
Subject: RE:  Converting Mathematics slides into PDF
Hi Bharat, 
> I read some of the posts on this topic, but I don't think 
> they address the problem I am having.
> 
> Is there a way to generate PDF file of Mathematica slides 
> such that what shows on the screen occupies one page or less 
> in the PDF file. For example, I have graphs which do not even 
> fill the screen completely, but do not fit in the page when I 
> make the PDF file (that is, only a part shows up on the page 
> and other part is cut out). Similarly, some of the text that 
> fits in half the screen fills up multiple pages in the PDF file.
I would approach this using the \Printout\ environment rather than
\SlideShow\.  
consisting of slides with different page lengths.  For what you want, this
would render some slides readable and others with very small text.  
With the Printout environment, you can set page breaks so that each slide \
is
readable.  You can then print your notebook to a PDF.
Dave.
===
Subject: Re: \No more memory available\ -- a recurring problem
>I took my MacBook Pro into the university computer support & repair
>not a hardware problem. I ran memtest (www.memtestosx.org/) on the new
>memory stick alone (it's from OMNI Technology) and found no errors.
>That leaves me with three possibilities:
>1) my code can simply be too difficult for a personal computer and some
>sort of grid or supercomputer is required 
>2) there is a bug in Mathematica, in the act of porting Mathematica to
>Intel (i.e. to being universal), or in OSX which causes memory failures
>even when there is plenty of memory and virtual memory available 
>3) my code is poorly written and needs radical change. 
You might try posting some of the code here. With the code, someone might be \
able to offer suggestions as to a work around, or at the very least give you \
an indication of which of the three possibilities above is most likely. Of \
course if your problem is so complex that it literally requires a few 100 \
lines of code, it is most likely not going to get much response nor is it \
really feasible to post such large code listings here.
--
To reply via email subtract one hundred and four
===
Subject: Re: Something in the spirit of letrec
>I'm trying to somehow come up with a function WithSelf such that
>evaluating
>WithSelf[{a,self[1]+b,self[2]+2}]
>will result in {a,a+b,a+b+c}. 
FoldList does this, i.e.,
In[2]:=
Rest@FoldList[Plus,0,{a,b,c}]
Out[2]=
{a,a+b,a+b+c}
--
To reply via email subtract one hundred and four
===
Subject: Re: Re:  \No more memory available\ -- a recurring problem
Eureka! I figured it out. In short, I had the NDSolve MaxSteps option
set to Infinity, which was causing the \No more memory available\
errors. By reducing this option to 20,000 and wrapping Check[] around
the instances in which I use NDSolve, I seem to have solved the
problem.
-- J\.87nos Lobb, Bill Rowe, James Gilmore, and Bruce Miller. Your
thoughtful comments and persistence uplifted me from some dark days. I
can't thank you enough.
Now, finally, back to the research...
Charlie
> Charlie,
 I would set $RecirsionLimit and its cousins to Infinity.  You might
> also hit a built in limit of Mathematica that is not obvious and
> might be processor dependent.  /I managed to to that with GraphPlot
> and the G5 duo Core :) /
 Are you able to run other Mathematica programs on the MacBook Pro
> fine ?  /Like creating a huge array that would eat up minimum half a
> GB memory/.  If yes, then you have to look your program with a fine
> tooth comb and nail it down where does it fail.  You might able to do
> it with a group of Print statements placed at the right places in
> your program.  If you can nail down what particular statement is
> offending the MacBook Pro, then you will have a good chance to fix it.
 If it dies in the middle of an NDsolve, then try that statement with
> different parameters if you can or try to program it out with a
> replacement/approximation.
 You might want to look the IMTEK Mathematica Supplement (IMS)
 http://www.imtek.de/simulation/
 where they provide some other methods to solve PDE-s and use their
> packages to see if those fail or not.
 Good luck,
 J\.87nos


 J\.87nos,
 Below are my responses.
> Do you have any Mathematica.crash.log or Mathkernel.crash.log file in
>> your ~Library/Logs/Crashreporter folder ?
>> If yes, does the datetime entries coincide with the kernel crashes ?
 I don't see any Mathematica crash logs in my
> ~Library/Logs/Crashreporter folder (only one file called
> mds.crash.log, which I doubt is Mathematica-related).
> Did you try to eliminate or cut back on History /$HistoryLength /?
 I had set $HistoryLength to 5 all along. Since learning that
> $HistoryLength can take up a lot of memory, I set it to zero, but I
> still get the \No more memory available\ message often.
> Did you set any iteration related env variables Wolfram provide to
>> Infinity like $IterationLimit and $RecusionLimit, etc....
 No, I did not change any of those from the default. I rarely, if ever,
> get other error messages, such as \reached the recursion limit,\
> preceding the \No more memory available\ message.

> Charlie

>> I took my MacBook Pro into the university computer support & repair
>> probably
>> is not a hardware problem. I ran memtest (www.memtestosx.org/)
>> on the
>> new memory stick alone (it's from OMNI Technology) and found no
>> errors.
>> That leaves me with three possibilities:
>> 1) my code can simply be too difficult for a personal computer and
>> some sort of grid or supercomputer is required
>> 2) there is a bug in Mathematica, in the act of porting
>> Mathematica to
>> Intel (i.e. to being universal), or in OSX which causes memory
>> failures even when there is plenty of memory and virtual memory
>> available
>> 3) my code is poorly written and needs radical change. To answer
>> Bill
>> Rowe's question, I would consider myself an intermediate
>> Mathematica
>> user, and I think I'm fairly well-versed in the area of
>> Mathematica in
>> which I'm working (namely, solving and analyzing differential
>> equations).
>> Basically, my calculation consists of using NDSolve thousands of
>> times
>> to solve and evaluate PDE's in order to find ones that are chaotic.
>> I'm not sure what else to do. Contact Wolfram maybe?
>> I guess this ends my research project.
>> Charlie
> Charlie,
 Your MacBook Pro has an Intel  core duo processor and looks to me
> that Mathematica is using just one of the cores.  I do not know
>> what
> You can do about it.  Probably nothing.  If you can put your
>> hand on
> a Personal Grid edition that would engage the second core more
> vigorously.
 I do not see that you have a memory problem as regarding the
>> usage of
> the memory reported by top.
 You might  have a bad memory chip and that is crashing your
> program.   Unfortunately we see bad memory coming even from Apple
> these days with newer Macs.  There is a service CD that came with
> your MacBook Pro.  Try to run the memory test from it all night
>> in a
> loop and see if it is caching anything.
 You might also want to do a hard drive check with Diskwarrior.
>> Wooly
> disks can lead to wooly virtual memory and causing crashes.
 Look if you have anything in ~/Library/Logs/Crashreporter
>> regarding
> Mathematica.crash.log or Mathkernel.cash.log.   /I have plenty in
> mine for both :)/  If you find any, send them to
>> support@wolfram.com
> and ask their opinion.
 I also recommend to do additional exhausting tests for the
> motherboard and other elements of the machine by using the service
> CD.
 In summary it looks like a hardware problem - unless you have
>> other
> long running Mathematica programs which run just fine.
 With the best,
 J\.87nos

>> J\.87nos,
>> information in the terminal window, so I have attached
>> screenshots of
>> it at various stages (see the names of the images). I have
>> attached
>> them in chronological order.
>> It seems to me that Mathematica is only being allocated 50% of
>> the CPU
>> PowerBook,
>> 93% of the CPU was being used while Mathematica was running my
>> code.
>> How do I change this on my 15\ MacBook Pro?
>> My 12\ PowerBook is able to run this code overnight and very
>> rarely
>> gives the \ran out of memory\ error message. I would like to
>> be able
>> to run my code overnight on my MacBook Pro, but it can only run
>> for 5
>> to 10 minutes before it runs out of memory and completely halts.
>> Thus
>> I can't \set it and leave it\; instead, I must always be at my
>> computer to restart the calculation every 5 to 10 minutes. Any
>> suggestions?
>> Charlie
>>>> I keep getting the following error message while running a
>> search
>> of many PDE's:
>> \No more memory available.
>> Mathematica kernel has shut down.
>> Try quitting other applications and then retry.\
>> This message is becoming so common that it is crippling my
>> research
>> project, which is to find the simplest PDE with one quadratic
>> nonlinearity that is chaotic.
>> I have Googled and searched the Mathgroup archives for help,
>> and I
>> employed the following fixes:
>> --CODE & SOFTWARE--
>> 1) Share[] -- does not help because my code rarely has common
>>> elements
>> that Share[] could consolidate.
>> 2) I used Module[] and made as many variables local as
>> possible. I
>> eliminated extra variables and functions in addition to adding
>>> Clear[]
>> in several places in the code to clear variables that are no
>> longer
>> needed.
>> 3) I streamlined and optimized the code in general, and I
>> made the
>> routines as simple and least data-intensive as possible.
>> 4) I never have any other applications open when running
>>> Mathematica.
>> (Unfortunately, closing other applications is the only
>>> suggestion that
>> Mathematica provides in the \out of memory\ error message.)
>> 5) I should note that I run the latest version of Mathematica
>> (5.2.2.0).
>> --HARDWARE--
>> 1) I upgraded from a 12\ PowerBook to a 15\ MacBook Pro. This
>>> computer
>> has an 2.16 GHz Intel Core Duo processor.
>> 3) I now have much more free hard drive space (20 GB) in case
>> Mathematica needs to use virtual memory. In addition, this hard
>>> drive
>> is 7200 rpm versus the standard 5400 rpm.
>> I have been using functions like MemoryInUse[] and On
>>> [MemoryConserve]
>> to monitor the use of memory. I almost always find that I am
>> only
>> using a few megabytes at a time (usually 5-10 MB, sometimes as
>>> high as
>> 100 or once 1000 after a large computation). However, despite
>>> finding
>> that I usually use only a few megabytes during my
>> computations, I
>> often get the above \out of memory\ error message, and I can
>> never
>> figure out why and how much memory was needed in that
>> particular
>> computation.
>> Some questions:
>> 1) How can I tell whether Mathematica is using virtual
>> memory or
>>> not?
>> It seems to me that it is not. I have 20 GB of hard drive space
>>> free
>> and it never seems to use it. Why should I ever run out of
>>> memory if I
>> provides?
>>> Mathematica can
>> access?
>>> greatly
>> appreciate it. If you need more information on my code or the
>> computations I am doing, I would be happy to provide it.
>> Charlie
>>>> If you open a Terminal window on the side and type \top\ and \
hit
>>> return and start up your application, what do You see for real
>>> memory
>>> usage for the kernel and the front end from top?
>>>> J\.87nos
>>> <before mathematica.jpeg>> <mma loaded, code loaded, pre-execution.jpeg>> <running1.jpeg>> <running2.jpeg>> <running3.jpeg>> <running4.jpeg>> <ran out of memory.jpeg


> ----------------------------------------------
> Trying to argue with a politician is like lifting up the head of a
> corpse.
> (S. Lem: His Master Voice)

===
Subject: Re: Re:  \No more memory available\ -- a recurring problem
J\.87nos,
Below are my responses.
> Do you have any Mathematica.crash.log or Mathkernel.crash.log file in
> your ~Library/Logs/Crashreporter folder ?
> If yes, does the datetime entries coincide with the kernel crashes ?
I don't see any Mathematica crash logs in my
~Library/Logs/Crashreporter folder (only one file called
mds.crash.log, which I doubt is Mathematica-related).
> Did you try to eliminate or cut back on History /$HistoryLength /?
I had set $HistoryLength to 5 all along. Since learning that
$HistoryLength can take up a lot of memory, I set it to zero, but I
still get the \No more memory available\ message often.
> Did you set any iteration related env variables Wolfram provide to
> Infinity like $IterationLimit and $RecusionLimit, etc....
No, I did not change any of those from the default. I rarely, if ever,
get other error messages, such as \reached the recursion limit,\
preceding the \No more memory available\ message.
Charlie
 I took my MacBook Pro into the university computer support & repair
> is not a hardware problem. I ran memtest (www.memtestosx.org/) on the
> new memory stick alone (it's from OMNI Technology) and found no
> errors.
 That leaves me with three possibilities:
> 1) my code can simply be too difficult for a personal computer and
> some sort of grid or supercomputer is required
> 2) there is a bug in Mathematica, in the act of porting Mathematica to
> Intel (i.e. to being universal), or in OSX which causes memory
> failures even when there is plenty of memory and virtual memory
> available
> 3) my code is poorly written and needs radical change. To answer Bill
> Rowe's question, I would consider myself an intermediate Mathematica
> user, and I think I'm fairly well-versed in the area of Mathematica in
> which I'm working (namely, solving and analyzing differential
> equations).
 Basically, my calculation consists of using NDSolve thousands of times
> to solve and evaluate PDE's in order to find ones that are chaotic.
 I'm not sure what else to do. Contact Wolfram maybe?
 I guess this ends my research project.

> Charlie
> Charlie,
>> Your MacBook Pro has an Intel  core duo processor and looks to me
>> that Mathematica is using just one of the cores.  I do not know what
>> You can do about it.  Probably nothing.  If you can put your hand on
>> a Personal Grid edition that would engage the second core more
>> vigorously.
>> I do not see that you have a memory problem as regarding the usage of
>> the memory reported by top.
>> You might  have a bad memory chip and that is crashing your
>> program.   Unfortunately we see bad memory coming even from Apple
>> these days with newer Macs.  There is a service CD that came with
>> your MacBook Pro.  Try to run the memory test from it all night in a
>> loop and see if it is caching anything.
>> You might also want to do a hard drive check with Diskwarrior.  Wooly
>> disks can lead to wooly virtual memory and causing crashes.
>> Look if you have anything in ~/Library/Logs/Crashreporter regarding
>> Mathematica.crash.log or Mathkernel.cash.log.   /I have plenty in
>> mine for both :)/  If you find any, send them to support@wolfram.com
>> and ask their opinion.
>> I also recommend to do additional exhausting tests for the
>> motherboard and other elements of the machine by using the service
>> CD.
>> In summary it looks like a hardware problem - unless you have other
>> long running Mathematica programs which run just fine.
>> With the best,
>> J\.87nos
> J\.87nos,
 information in the terminal window, so I have attached
> screenshots of
> it at various stages (see the names of the images). I have attached
> them in chronological order.
 It seems to me that Mathematica is only being allocated 50% of
> the CPU
> 93% of the CPU was being used while Mathematica was running my code.
> How do I change this on my 15\ MacBook Pro?
 My 12\ PowerBook is able to run this code overnight and very rarely
> gives the \ran out of memory\ error message. I would like to be \
able
> to run my code overnight on my MacBook Pro, but it can only run
> for 5
> to 10 minutes before it runs out of memory and completely halts.
> Thus
> I can't \set it and leave it\; instead, I must always be at my
> computer to restart the calculation every 5 to 10 minutes. Any
> suggestions?

 Charlie
>> I keep getting the following error message while running a search
>>> of many PDE's:
>>>> \No more memory available.
>>> Mathematica kernel has shut down.
>>> Try quitting other applications and then retry.\
>>>> This message is becoming so common that it is crippling my
>>> research
>>> project, which is to find the simplest PDE with one quadratic
>>> nonlinearity that is chaotic.
>>>> I have Googled and searched the Mathgroup archives for help, and I
>>> employed the following fixes:
>>>> --CODE & SOFTWARE--
>>> 1) Share[] -- does not help because my code rarely has common
>> elements
>>> that Share[] could consolidate.
>>>> 2) I used Module[] and made as many variables local as possible. I
>>> eliminated extra variables and functions in addition to adding
>> Clear[]
>>> in several places in the code to clear variables that are no
>>> longer
>>> needed.
>>>> 3) I streamlined and optimized the code in general, and I made the
>>> routines as simple and least data-intensive as possible.
>>>> 4) I never have any other applications open when running
>> Mathematica.
>>> (Unfortunately, closing other applications is the only
>> suggestion that
>>> Mathematica provides in the \out of memory\ error message.)
>>>> 5) I should note that I run the latest version of Mathematica
>>> (5.2.2.0).
>>>> --HARDWARE--
>>> 1) I upgraded from a 12\ PowerBook to a 15\ MacBook Pro. This
>> computer
>>> has an 2.16 GHz Intel Core Duo processor.
>>>>> 3) I now have much more free hard drive space (20 GB) in case
>>> Mathematica needs to use virtual memory. In addition, this hard
>> drive
>>> is 7200 rpm versus the standard 5400 rpm.
>>>>> I have been using functions like MemoryInUse[] and On
>> [MemoryConserve]
>>> to monitor the use of memory. I almost always find that I am only
>>> using a few megabytes at a time (usually 5-10 MB, sometimes as
>> high as
>>> 100 or once 1000 after a large computation). However, despite
>> finding
>>> that I usually use only a few megabytes during my computations, I
>>> often get the above \out of memory\ error message, and I can \
never
>>> figure out why and how much memory was needed in that particular
>>> computation.
>>>>> Some questions:
>>> 1) How can I tell whether Mathematica is using virtual memory or
>> not?
>>> It seems to me that it is not. I have 20 GB of hard drive space
>> free
>>> and it never seems to use it. Why should I ever run out of
>> memory if I
>>> Mathematica can
>>> access?
>>> greatly
>>> appreciate it. If you need more information on my code or the
>>> computations I am doing, I would be happy to provide it.
>>>> Charlie
>> If you open a Terminal window on the side and type \top\ and hit
>> return and start up your application, what do You see for real
>> memory
>> usage for the kernel and the front end from top?
>> J\.87nos
>> <before mathematica.jpeg> <mma loaded, code loaded, pre-execution.jpeg> <running1.jpeg> <running2.jpeg> <running3.jpeg> <running4.jpeg> <ran out of memory.jpeg
>
===
Subject: Re:  \No more memory available\ -- a recurring problem
Charlie
Do you have any Mathematica.crash.log or Mathkernel.crash.log file in  
your ~Library/Logs/Crashreporter folder ?
If yes, does the datetime entries coincide with the kernel crashes ?
If you find any of the two files, send them to support@wolfram.com.
Did you try to eliminate or cut back on History /$HistoryLength /?
Did you set any iteration related env variables Wolfram provide to  
Infinity like $IterationLimit and $RecusionLimit, etc....
With the best,
J\.87nos
> I took my MacBook Pro into the university computer support & repair
> is not a hardware problem. I ran memtest (www.memtestosx.org/) on the
> new memory stick alone (it's from OMNI Technology) and found no
> errors.
 That leaves me with three possibilities:
> 1) my code can simply be too difficult for a personal computer and
> some sort of grid or supercomputer is required
> 2) there is a bug in Mathematica, in the act of porting Mathematica to
> Intel (i.e. to being universal), or in OSX which causes memory
> failures even when there is plenty of memory and virtual memory
> available
> 3) my code is poorly written and needs radical change. To answer Bill
> Rowe's question, I would consider myself an intermediate Mathematica
> user, and I think I'm fairly well-versed in the area of Mathematica in
> which I'm working (namely, solving and analyzing differential
> equations).
 Basically, my calculation consists of using NDSolve thousands of times
> to solve and evaluate PDE's in order to find ones that are chaotic.
 I'm not sure what else to do. Contact Wolfram maybe?
 I guess this ends my research project.

> Charlie
> Charlie,
>> Your MacBook Pro has an Intel  core duo processor and looks to me
>> that Mathematica is using just one of the cores.  I do not know what
>> You can do about it.  Probably nothing.  If you can put your hand on
>> a Personal Grid edition that would engage the second core more
>> vigorously.
>> I do not see that you have a memory problem as regarding the usage of
>> the memory reported by top.
>> You might  have a bad memory chip and that is crashing your
>> program.   Unfortunately we see bad memory coming even from Apple
>> these days with newer Macs.  There is a service CD that came with
>> your MacBook Pro.  Try to run the memory test from it all night in a
>> loop and see if it is caching anything.
>> You might also want to do a hard drive check with Diskwarrior.  Wooly
>> disks can lead to wooly virtual memory and causing crashes.
>> Look if you have anything in ~/Library/Logs/Crashreporter regarding
>> Mathematica.crash.log or Mathkernel.cash.log.   /I have plenty in
>> mine for both :)/  If you find any, send them to support@wolfram.com
>> and ask their opinion.
>> I also recommend to do additional exhausting tests for the
>> motherboard and other elements of the machine by using the service  
>> CD.
>> In summary it looks like a hardware problem - unless you have other
>> long running Mathematica programs which run just fine.
>> With the best,
>> J\.87nos
> J\.87nos,
 information in the terminal window, so I have attached  
> screenshots of
> it at various stages (see the names of the images). I have attached
> them in chronological order.
 It seems to me that Mathematica is only being allocated 50% of  
> the CPU
> 93% of the CPU was being used while Mathematica was running my code.
> How do I change this on my 15\ MacBook Pro?
 My 12\ PowerBook is able to run this code overnight and very rarely
> gives the \ran out of memory\ error message. I would like to be able
> to run my code overnight on my MacBook Pro, but it can only run  
> for 5
> to 10 minutes before it runs out of memory and completely halts.  
> Thus
> I can't \set it and leave it\; instead, I must always be at my
> computer to restart the calculation every 5 to 10 minutes. Any
> suggestions?

 Charlie
>> I keep getting the following error message while running a search
>>> of many PDE's:
>>>> \No more memory available.
>>> Mathematica kernel has shut down.
>>> Try quitting other applications and then retry.\
>>>> This message is becoming so common that it is crippling my  
>>> research
>>> project, which is to find the simplest PDE with one quadratic
>>> nonlinearity that is chaotic.
>>>> I have Googled and searched the Mathgroup archives for help, and I
>>> employed the following fixes:
>>>> --CODE & SOFTWARE--
>>> 1) Share[] -- does not help because my code rarely has common
>> elements
>>> that Share[] could consolidate.
>>>> 2) I used Module[] and made as many variables local as possible. I
>>> eliminated extra variables and functions in addition to adding
>> Clear[]
>>> in several places in the code to clear variables that are no  
>>> longer
>>> needed.
>>>> 3) I streamlined and optimized the code in general, and I made the
>>> routines as simple and least data-intensive as possible.
>>>> 4) I never have any other applications open when running
>> Mathematica.
>>> (Unfortunately, closing other applications is the only
>> suggestion that
>>> Mathematica provides in the \out of memory\ error message.)
>>>> 5) I should note that I run the latest version of Mathematica
>>> (5.2.2.0).
>>>> --HARDWARE--
>>> 1) I upgraded from a 12\ PowerBook to a 15\ MacBook Pro. This
>> computer
>>> has an 2.16 GHz Intel Core Duo processor.
>>>>> 3) I now have much more free hard drive space (20 GB) in case
>>> Mathematica needs to use virtual memory. In addition, this hard
>> drive
>>> is 7200 rpm versus the standard 5400 rpm.
>>>>> I have been using functions like MemoryInUse[] and On
>> [MemoryConserve]
>>> to monitor the use of memory. I almost always find that I am only
>>> using a few megabytes at a time (usually 5-10 MB, sometimes as
>> high as
>>> 100 or once 1000 after a large computation). However, despite
>> finding
>>> that I usually use only a few megabytes during my computations, I
>>> often get the above \out of memory\ error message, and I can never
>>> figure out why and how much memory was needed in that particular
>>> computation.
>>>>> Some questions:
>>> 1) How can I tell whether Mathematica is using virtual memory or
>> not?
>>> It seems to me that it is not. I have 20 GB of hard drive space
>> free
>>> and it never seems to use it. Why should I ever run out of
>> memory if I
>>> Mathematica can
>>> access?
>>> greatly
>>> appreciate it. If you need more information on my code or the
>>> computations I am doing, I would be happy to provide it.
>>>> Charlie
>> If you open a Terminal window on the side and type \top\ and hit
>> return and start up your application, what do You see for real  
>> memory
>> usage for the kernel and the front end from top?
>> J\.87nos
>> <before mathematica.jpeg> <mma loaded, code loaded, pre-execution.jpeg> <running1.jpeg> <running2.jpeg> <running3.jpeg> <running4.jpeg> <ran out of memory.jpeg>
===
Subject: Re:  \No more memory available\ -- a recurring problem
I took my MacBook Pro into the university computer support & repair
is not a hardware problem. I ran memtest (www.memtestosx.org/) on the
new memory stick alone (it's from OMNI Technology) and found no
errors.
That leaves me with three possibilities:
1) my code can simply be too difficult for a personal computer and
some sort of grid or supercomputer is required
2) there is a bug in Mathematica, in the act of porting Mathematica to
Intel (i.e. to being universal), or in OSX which causes memory
failures even when there is plenty of memory and virtual memory
available
3) my code is poorly written and needs radical change. To answer Bill
Rowe's question, I would consider myself an intermediate Mathematica
user, and I think I'm fairly well-versed in the area of Mathematica in
which I'm working (namely, solving and analyzing differential
equations).
Basically, my calculation consists of using NDSolve thousands of times
to solve and evaluate PDE's in order to find ones that are chaotic.
I'm not sure what else to do. Contact Wolfram maybe?
I guess this ends my research project.
Charlie
> Charlie,
 Your MacBook Pro has an Intel  core duo processor and looks to me
> that Mathematica is using just one of the cores.  I do not know what
> You can do about it.  Probably nothing.  If you can put your hand on
> a Personal Grid edition that would engage the second core more
> vigorously.
 I do not see that you have a memory problem as regarding the usage of
> the memory reported by top.
 You might  have a bad memory chip and that is crashing your
> program.   Unfortunately we see bad memory coming even from Apple
> these days with newer Macs.  There is a service CD that came with
> your MacBook Pro.  Try to run the memory test from it all night in a
> loop and see if it is caching anything.
 You might also want to do a hard drive check with Diskwarrior.  Wooly
> disks can lead to wooly virtual memory and causing crashes.
 Look if you have anything in ~/Library/Logs/Crashreporter regarding
> Mathematica.crash.log or Mathkernel.cash.log.   /I have plenty in
> mine for both :)/  If you find any, send them to support@wolfram.com
> and ask their opinion.
 I also recommend to do additional exhausting tests for the
> motherboard and other elements of the machine by using the service CD.
 In summary it looks like a hardware problem - unless you have other
> long running Mathematica programs which run just fine.
 With the best,
 J\.87nos

> J\.87nos,
 information in the terminal window, so I have attached screenshots of
> it at various stages (see the names of the images). I have attached
> them in chronological order.
 It seems to me that Mathematica is only being allocated 50% of the CPU
> 93% of the CPU was being used while Mathematica was running my code.
> How do I change this on my 15\ MacBook Pro?
 My 12\ PowerBook is able to run this code overnight and very rarely
> gives the \ran out of memory\ error message. I would like to be able
> to run my code overnight on my MacBook Pro, but it can only run for 5
> to 10 minutes before it runs out of memory and completely halts. Thus
> I can't \set it and leave it\; instead, I must always be at my
> computer to restart the calculation every 5 to 10 minutes. Any
> suggestions?

 Charlie
> I keep getting the following error message while running a search
>> of many PDE's:
>> \No more memory available.
>> Mathematica kernel has shut down.
>> Try quitting other applications and then retry.\
>> This message is becoming so common that it is crippling my research
>> project, which is to find the simplest PDE with one quadratic
>> nonlinearity that is chaotic.
>> I have Googled and searched the Mathgroup archives for help, and I
>> employed the following fixes:
>> --CODE & SOFTWARE--
>> 1) Share[] -- does not help because my code rarely has common
>> elements
>> that Share[] could consolidate.
>> 2) I used Module[] and made as many variables local as possible. I
>> eliminated extra variables and functions in addition to adding
>> Clear[]
>> in several places in the code to clear variables that are no longer
>> needed.
>> 3) I streamlined and optimized the code in general, and I made the
>> routines as simple and least data-intensive as possible.
>> 4) I never have any other applications open when running
>> Mathematica.
>> (Unfortunately, closing other applications is the only
>> suggestion that
>> Mathematica provides in the \out of memory\ error message.)
>> 5) I should note that I run the latest version of Mathematica
>> (5.2.2.0).
>> --HARDWARE--
>> 1) I upgraded from a 12\ PowerBook to a 15\ MacBook Pro. This
>> computer
>> has an 2.16 GHz Intel Core Duo processor.
>> 3) I now have much more free hard drive space (20 GB) in case
>> Mathematica needs to use virtual memory. In addition, this hard
>> drive
>> is 7200 rpm versus the standard 5400 rpm.
>> I have been using functions like MemoryInUse[] and On
>> [MemoryConserve]
>> to monitor the use of memory. I almost always find that I am only
>> using a few megabytes at a time (usually 5-10 MB, sometimes as
>> high as
>> 100 or once 1000 after a large computation). However, despite
>> finding
>> that I usually use only a few megabytes during my computations, I
>> often get the above \out of memory\ error message, and I can never
>> figure out why and how much memory was needed in that particular
>> computation.
>> Some questions:
>> 1) How can I tell whether Mathematica is using virtual memory or
>> not?
>> It seems to me that it is not. I have 20 GB of hard drive space
>> free
>> and it never seems to use it. Why should I ever run out of
>> memory if I
>> Mathematica can
>> access?
>> greatly
>> appreciate it. If you need more information on my code or the
>> computations I am doing, I would be happy to provide it.
>> Charlie
>> If you open a Terminal window on the side and type \top\ and hit
>> return and start up your application, what do You see for real memory
>> usage for the kernel and the front end from top?
>> J\.87nos
>> <before mathematica.jpeg> <mma loaded, code loaded, pre-execution.jpeg> <running1.jpeg> <running2.jpeg> <running3.jpeg> <running4.jpeg> <ran out of memory.jpeg
>
===
Subject: Re:  Re: Finding the Number of Pythagorean Triples below a bound
becomes really fast ;-)
countT = Compile[{{m, _Integer}, {k, _Integer}}, Module[
   {i = 0, c2, diff, sdiff},
   Do [
    If[Mod[c^2 + k, 4] != 3, c2 = Floor[Sqrt[(c^2 + k)/2]];
     Do [
       diff = c^2 + k - b^2;
          sdiff = Sqrt[N[diff]];
          If [sdiff >= b && sdiff == Round[sdiff], i++];
       , {b, 1, c2}]];
, {c, 1, 10^m - 1}];
i]]
countT[3,2]+countT[3,-2]//Timing
{1.0032 Second,282}
I will try the case m=4 and m=5 and send you the results. I am not  
promising to do this very soon, just in case ;-)
Andrzej Kozlowski
> Here is the fastest code I have seen so far that works (it is based  
> on the one that Daniel sent you). I have corrected and enhanced   
> and enhanced it. It gives approximate answers for reasons that I  
> explained in earlier postings (use of machine arithmetic to test  
> for perfect integers). Of course it is easy to replace the code by  
> exact code by replacing the numerical test for a perfect square by  
> an exact one.

> countTriples[m_,k_] := Module[
>   {i=0, c2, diff, sdiff},
>   Do [
>    If[Mod[c^2+k,4]!=3, c2 = Floor[Sqrt[(c^2+k)/2]];
>     Do [
>       diff = c^2+k-b^2;
>          sdiff = Sqrt[N[diff]];
>          If [sdiff>=b&&sdiff==Round[sdiff],i++];
>       , {b,1,c2}]];
> ,{c,1,10^m-1}];
> i]

> countTriples[3,2]+countTriples[3,-2]//Timing

> {12.3746 Second,282}
 The correct answer is 283.
 This code should easily deal with the case m=4 (I have not yet  
> tried it) and I think even m=5 should now be within reach.
 Andrzej Kozlowski


>> The \improvement\ below which I sent a little earlier is wrong  
>> (even though it returned correct answers). Obviously the point is  
>> that c^2+k  can't be of the form 4n + 3, but there is no reason  
>> why a^2+b^2-k can't be of that form. Since my code does not  
>> explicitly select c it can't make use of this additional  
>> improvement. A different code, which uses explicit choices of  
>> (say) a and c and tests for b being a perfect square could exploit  
>> this fact and perhaps gain extra speed. It should not be difficult  
>> to write such a code along the lines I have been using.
>> Andrzej
>
>>> Hello all,
>>>> emailed
>>> me.
>>>> To recall, the problem was counting the number of solutions to a
>>> bivariate polynomial equal to a square,
>>>> Poly(a,b) = c^2
>>>> One form that interested me was the Pythagorean-like equation:
>>>> a^2 + b^2 = c^2 + k
>>>> for {a,b} a positive integer, 0<a<=b, and k any small  
>>> integer. I was
>>> wondering about the density of solutions to this since I knew  
>>> in the
>>> special case of k=0, let S(N) be the number of primitive  
>>> solutions
>>> with
>>> c < N, then S(N)/N = 1/(2pi) as N -> inf.
>>>> For k a squarefree integer, it is convenient that any  
>>> solution is
>>> also
>>> primitive. I used a simple code that allowed me to find S 
>>> (10^m) with
>>> m=1,2,3 for small values of k (for m=4 took my code more than  
>>> 30 mins
>>> so I aborted it). The data is given below:
>>>> Note: Values are total S(N) for *both* k & -k:
>>>> k = 2
>>> S(N) = 4, 30, 283
>>>> k = 3
>>> S(N) = 3, 41, 410
>>>> k = 5
>>> S(N) = 3, 43, 426
>>>> k = 6
>>> S(N) = 3, 36, 351
>>>> Question: Does S(N)/N for these also converge? For example,  
>>> for the
>>> particular case of k = -6, we have
>>>> S(N) = 2, 20, 202
>>>> which looks suspiciously like the ratio might be converging.
>>>> Anybody know of a code for this that can find m=4,5,6 in a  
>>> reasonable
>>> amount of time?
>>>>> Yours,
>>>> Titus
>>> Here is a piece code which utilises the ideas I have described in
>> my previous posts:
>> ls = Prime /@ Range[3, 10];
>> test[n_] :=
>>       Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>> Integers]
>> f[P_, k_] := Sum[If[(w =
>>    a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a,
>>           Floor[Sqrt[P^2 - a^2]]}]
>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>> We can easily confirm the results of your computations,.e.g.  
>> for k=2.
>> Table[g[i,2],{i,3}]
>> {4,30,283}
>> Since you have not revealed the \simple code\ you have used it is
>> hard to tell if the above one is any better. It is however,
>> certainly capable of solving the problem for m=4:
>> g[4,2]//Timing
>> {4779.39 Second,2763}
>> The time it took on my 1 Gigahertz PowerBook was over 70 minutes,
>> which is longer than you thought \reasonable\, so I am still not
>> sure if this is any improvement on what you already have. The  
>> time
>> complexity of this algorithm seems somewhat larger than  
>> exponential
>> so I would expect that it will take about 6 hours to deal with  
>> n=5
>> on my computer, and perhaps 2 weeks to deal with n=6.
>> Andrzej Kozlowski
>>>>> I mistakenly copied and pasted a wrong (earlier) definition of f.
>>> Here is the correct one:
>>>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>      1 && (w = a^2 +
>>>           b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>>            P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]]
>>>> The definition of g is as before:
>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>> Andrzej Kozlowski
>>> Below is a faster version of the above code. (It owes a  
>> significant improvement to Daniel Lichtblau, which I borrowed  
>> from him without his knowledge ;-))
>> test[n_] :=
>>     JacobiSymbol[n,
>>   Prime[Random[Integer, {2, 20}]]] \[CapitalEth] -1 && \
Element[Sqrt[n],
>> Integers]
>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>     1 && (w = a^2 +
>>          b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>> +k]]}]]
>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>> The improvement is in the upper bound on a in the sum. Since a  
>> is the smaller of the two squares whose sum is equal to P^2+k it  
>> can't be larger than  Floor[Sqrt[(P^2+k)/2]].
>> Note that you can improve the performance by loosing some  
>> accuracy if you use a cruder test for a perfect square:
>> test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w]
>>    ]
>> f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>     1 && (w = a^2 +
>>          b^2 - k) < P^2 && test1[w], 1, 0], {a, 1,
>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>> +k]]}]]
>> g1[m_, k_] := f1[10^m, k] + f1[10^m, -k]
>> Let's compare the two cases.
>> In[7]:=
>> g[3,2]//Timing
>> Out[7]=
>> {89.554 Second,283}
>> In[8]:=
>> g1[3,2]//Timing
>> Out[8]=
>> {37.376 Second,283}
>> So we see that we get the same answer and the second approach is  
>> considerably faster. However:
>> In[9]:=
>> g[1,6]//Timing
>> Out[9]=
>> {0.008863 Second,3}
>> In[10]:=
>> g1[1,6]//Timing
>> Out[10]=
>> {0.005429 Second,5}
>> The correct answer is 3 (returned by the first method). The  
>> faster method found two false solutions. This should not matter  
>> if you are interested only in approximate answers (as you seem  
>> to be) but it is worth keeping in mind.
>> Andrzej Kozlowski

> I have noticed one more obvious improvement. We can replace test by:
 test[n_] :=
>     Mod[n, 4] =!= 3 && JacobiSymbol[n,
>   Prime[Random[Integer, {2, 100}]]] \[CapitalEth] -1 && \
Element[Sqrt[n],
> Integers]
 and test1 by
 test1[n_] := With[{w = Sqrt[N[n]]}, Mod[n, 4] =!= 3 && w == Round[w]
>    ]
 We are simply making use of the easily to prove fact that an  
> integer of the form 4 k + 1 can never be the sum of two squares.  
> There is a noticeable improvement in the performance of g:

> g[3,2]//Timing
 {58.0786 Second,283}
 However, the performance of g1 seems to actually decline slightly:

> g1[3,2]//Timing
 {40.8776 Second,283}

> However, there are fewer cases of \false solutions\:
 In[22]:=
> g1[1,6]//Timing
 Out[22]=
> {0.006694 Second,4}

> I am still not sure what is the most efficient use of  
> JacobiSymbol in this kind of problems. In the first version of my  
> code I used a test involving the first few odd primes, afterwards  
> I switched to using just one random prime in taken form a certain  
> range. Evaluation of JacobiSympol[m,p] != -1 is certainly much  
> faster than that of Element[Sqrt[m],Integers], but it is not  
> clear what is the most efficient number of primes to use, which  
> are the best primes and whether it is better to choose them at  
> random or just use a fixed selection.
> The numerical test Sqrt[N[n]]==Round[Sqrt[N[n]] is even faster,  
> but it will sometimes produce \false squares\.
 Andrzej Kozlowski

>
===
Subject: Re:  Re: Finding the Number of Pythagorean Triples below a bound
Daniel Lichtblau has already informed me that running the code below  
on a 64 bit computer will not help at all, so there is no point  
trying :-(
However, he also suggested a way around the Compile problem with  
integers larger than machine integers :-)
I am trying it out right now.
Andrzej Kozlowski
> I need to add some corrections. There were some some small mistakes  
> in the code I posted earlier, which caused problems with running  
> the compiled code for negative values of k, and which also probably  
> accounted for the slightly incorrect answers which the  
> \approximate\ code returned. I attributed the inaccuracy to the use  
> of numerical precision in the test for a number being a perfect  
> square but it seems that the cause was (probably) elsewhere. I  
> don't want to try to make this message too long so I won't bother  
> explaining what I think the mistakes were; but I will give what I  
> think is the correct code. I have decided to separate the code for  
> negative and positive values of k. The code for negative k works  
> also for positive k's but is slightly slower, due to the extra test  
> that needs to be performed. For this reason, and for the sake of  
> greater clarity I have decided to separate the two codes.
 The code for positive k:
 countTriplesP = Compile[{{m, _Integer}, {k, _Integer}}, Module[
>   {i = 0, c2, diff, sdiff},
>   Do [
>    If[Mod[c^2 + k, 4] != 3, c2 = Floor[Sqrt[(c^2 + k)/2]];
>     Do [
>       diff = c^2 + k - b^2;
>          sdiff = Sqrt[N[diff]];
>          If [sdiff >= b && sdiff == Round[sdiff], i++];
>       , {b, 1, c2}]];
> , {c, 1, 10^m - 1}];
> i]]
 The code for negative k:
 countTriplesN = Compile[{{m, _Integer}, {k, _Integer}}, Module[
>   {i = 1, c2, diff, sdiff},
>   Do [
>    If[Mod[c^2 + k, 4] != 3 && c^2 + k Ò 0, c2 = \
Floor[Sqrt[(c^2 +  
> k)/2]];
>     Do [
>       diff = c^2 + k - b^2;
>          sdiff = Sqrt[N[diff]];
>          If [sdiff >= b && sdiff == Round[sdiff], i++];
>       , {b, 1, c2}]];
> , {c, 1, 10^m - 1}];
> i]]
 Now we get:
 countTriplesP[1,6]+countTriplesN[1,-6]//Timing

> {0.000221 Second,3}

> countTriplesP[3,2]+countTriplesN[3,-2]//Timing

> {0.95186 Second,282}

> countTriplesP[4,2]+countTriplesN[4,-2]//Timing

> {95.2177 Second,2762}


> Note that these values are consistently less by one form the values  
> obtained by Titus and also by me using my earlier \exact\ code, but  
> actually I believe that this disparity was due to mistakes in the  
> earlier code. In any case, if we replace the \numerical\ test for a  
> perfect square by the much slower \exact\ test, the answers will be  
> the same, so the difference of 1 is certainly not due to the use of  
> numerical precision. Anyway, everything works fast and looks  
> perfectly satisfactory but then there is a snag. I decided to run  
> the code for m=5 and k=2, started it and went out for several  
> hours. When I came back I was rather disappointed to see that it  
> was still running and then I saw the message:
 countTriplesP[5,2]//Timing
 CompiledFunction::\cfn\ Numerical error  encountered at instruction  
> 10; proceeding with uncompiled evaluation.
 I assume the cause of this is that on 32 bit computers Compile  
> cannot deal with integers larger than 2^32 but we have:

> c=10^5;
 Log[2,c^2]//N
 33.2193
 I can't at the moment see any way around this problem except to run  
> uncompiled code (far too long) or try a 64 bit computer.  
> Unfortunately I do not have one and I don't think Titus has one, so  
> if any body has a 64-bit computer with a 64 bit version of  
> Mathematica installed, and a little spare time, I think both of us  
> would like to know how the above code performs in this setting. If  
> it works it will provide a nice bit of advertising for 64 bit  
> computing.
 Andrzej Kozlowski


> My latest code has just managed:
>> countT[4,2]+countT[4,-2]//Timing
>> {93.9638 Second,2762}
>> That is less than two minutes on a 1 gigahertz computer. The  
>> correct answer is actually 2763 (by my earlier computation using  
>> the exact test) so we have lost one solution but gained more than  
>> 50 fold improvement in performance!
>> The case m=5 is now certainly feasible, although I am not sure if  
>> I wish my not very powerful PowerBook to be occupied for so long,  
>> as I need to use Mathematica for other tasks. Perhaps I can now  
>> leave this to others.
>> Andrzej Kozlowski
> it becomes really fast ;-)
 countT = Compile[{{m, _Integer}, {k, _Integer}}, Module[
>   {i = 0, c2, diff, sdiff},
>   Do [
>    If[Mod[c^2 + k, 4] != 3, c2 = Floor[Sqrt[(c^2 + k)/2]];
>     Do [
>       diff = c^2 + k - b^2;
>          sdiff = Sqrt[N[diff]];
>          If [sdiff >= b && sdiff == Round[sdiff], i++];
>       , {b, 1, c2}]];
> , {c, 1, 10^m - 1}];
> i]]

> countT[3,2]+countT[3,-2]//Timing

> {1.0032 Second,282}
 I will try the case m=4 and m=5 and send you the results. I am  
> not promising to do this very soon, just in case ;-)
 Andrzej Kozlowski


>> Here is the fastest code I have seen so far that works (it is  
>> based on the one that Daniel sent you). I have corrected and  
>> enhanced  and enhanced it. It gives approximate answers for  
>> reasons that I explained in earlier postings (use of machine  
>> arithmetic to test for perfect integers). Of course it is easy  
>> to replace the code by exact code by replacing the numerical  
>> test for a perfect square by an exact one.
>> countTriples[m_,k_] := Module[
>>   {i=0, c2, diff, sdiff},
>>   Do [
>>    If[Mod[c^2+k,4]!=3, c2 = Floor[Sqrt[(c^2+k)/2]];
>>     Do [
>>       diff = c^2+k-b^2;
>>          sdiff = Sqrt[N[diff]];
>>          If [sdiff>=b&&sdiff==Round[sdiff],i++];
>>       , {b,1,c2}]];
>> ,{c,1,10^m-1}];
>> i]
>> countTriples[3,2]+countTriples[3,-2]//Timing
>> {12.3746 Second,282}
>> The correct answer is 283.
>> This code should easily deal with the case m=4 (I have not yet  
>> tried it) and I think even m=5 should now be within reach.
>> Andrzej Kozlowski
>>> The \improvement\ below which I sent a little earlier is wrong  
>>> (even though it returned correct answers). Obviously the point  
>>> is that c^2+k  can't be of the form 4n + 3, but there is no  
>>> reason why a^2+b^2-k can't be of that form. Since my code does  
>>> not explicitly select c it can't make use of this additional  
>>> improvement. A different code, which uses explicit choices of  
>>> (say) a and c and tests for b being a perfect square could  
>>> exploit this fact and perhaps gain extra speed. It should not  
>>> be difficult to write such a code along the lines I have been  
>>> using.
>>>> Andrzej
>>>>>>>>>>>>> Hello all,
>>>>>> privately
>>>> emailed
>>>> me.
>>>>>> To recall, the problem was counting the number of  
>>>> solutions to a
>>>> bivariate polynomial equal to a square,
>>>>>> Poly(a,b) = c^2
>>>>>> One form that interested me was the Pythagorean-like  
>>>> equation:
>>>>>> a^2 + b^2 = c^2 + k
>>>>>> for {a,b} a positive integer, 0<a<=b, and k any small  
>>>> integer. I was
>>>> wondering about the density of solutions to this since I  
>>>> knew in the
>>>> special case of k=0, let S(N) be the number of primitive  
>>>> solutions
>>>> with
>>>> c < N, then S(N)/N = 1/(2pi) as N -> inf.
>>>>>> For k a squarefree integer, it is convenient that any  
>>>> solution is
>>>> also
>>>> primitive. I used a simple code that allowed me to find S 
>>>> (10^m) with
>>>> m=1,2,3 for small values of k (for m=4 took my code more  
>>>> than 30 mins
>>>> so I aborted it). The data is given below:
>>>>>> Note: Values are total S(N) for *both* k & -k:
>>>>>> k = 2
>>>> S(N) = 4, 30, 283
>>>>>> k = 3
>>>> S(N) = 3, 41, 410
>>>>>> k = 5
>>>> S(N) = 3, 43, 426
>>>>>> k = 6
>>>> S(N) = 3, 36, 351
>>>>>> Question: Does S(N)/N for these also converge? For  
>>>> example, for the
>>>> particular case of k = -6, we have
>>>>>> S(N) = 2, 20, 202
>>>>>> which looks suspiciously like the ratio might be converging.
>>>>>> Anybody know of a code for this that can find m=4,5,6 in a  
>>>> reasonable
>>>> amount of time?
>>>>>>>> Yours,
>>>>>> Titus
>>>>>>> Here is a piece code which utilises the ideas I have  
>>> described in
>>> my previous posts:
>>>> ls = Prime /@ Range[3, 10];
>>>> test[n_] :=
>>>       Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt 
>>> [n],
>>> Integers]
>>>> f[P_, k_] := Sum[If[(w =
>>>    a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a,
>>>           Floor[Sqrt[P^2 - a^2]]}]
>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>> We can easily confirm the results of your  
>>> computations,.e.g. for k=2.
>>>>> Table[g[i,2],{i,3}]
>>>>> {4,30,283}
>>>>> Since you have not revealed the \simple code\ you have used  
>>> it is
>>> hard to tell if the above one is any better. It is however,
>>> certainly capable of solving the problem for m=4:
>>>>> g[4,2]//Timing
>>>>> {4779.39 Second,2763}
>>>> The time it took on my 1 Gigahertz PowerBook was over 70  
>>> minutes,
>>> which is longer than you thought \reasonable\, so I am  
>>> still not
>>> sure if this is any improvement on what you already have.  
>>> The time
>>> complexity of this algorithm seems somewhat larger than  
>>> exponential
>>> so I would expect that it will take about 6 hours to deal  
>>> with n=5
>>> on my computer, and perhaps 2 weeks to deal with n=6.
>>>> Andrzej Kozlowski
>>>>>>>>>>> I mistakenly copied and pasted a wrong (earlier) definition  
>>>> of f.
>>>> Here is the correct one:
>>>>>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>>      1 && (w = a^2 +
>>>>           b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>>>            P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]]
>>>>>> The definition of g is as before:
>>>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>>>> Andrzej Kozlowski
>>>>>>> Below is a faster version of the above code. (It owes a  
>>> significant improvement to Daniel Lichtblau, which I borrowed  
>>> from him without his knowledge ;-))
>>>> test[n_] :=
>>>     JacobiSymbol[n,
>>>   Prime[Random[Integer, {2, 20}]]] \[CapitalEth] -1 && \
Element[Sqrt[n],
>>> Integers]
>>>>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>     1 && (w = a^2 +
>>>          b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>>> +k]]}]]
>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>> The improvement is in the upper bound on a in the sum. Since  
>>> a is the smaller of the two squares whose sum is equal to P^2 
>>> +k it can't be larger than  Floor[Sqrt[(P^2+k)/2]].
>>>> Note that you can improve the performance by loosing some  
>>> accuracy if you use a cruder test for a perfect square:
>>>> test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w]
>>>    ]
>>>> f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>     1 && (w = a^2 +
>>>          b^2 - k) < P^2 && test1[w], 1, 0], {a, 1,
>>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>>> +k]]}]]
>>>>> g1[m_, k_] := f1[10^m, k] + f1[10^m, -k]
>>>>> Let's compare the two cases.
>>>> In[7]:=
>>> g[3,2]//Timing
>>>> Out[7]=
>>> {89.554 Second,283}
>>>> In[8]:=
>>> g1[3,2]//Timing
>>>> Out[8]=
>>> {37.376 Second,283}
>>>> So we see that we get the same answer and the second approach  
>>> is considerably faster. However:
>>>> In[9]:=
>>> g[1,6]//Timing
>>>> Out[9]=
>>> {0.008863 Second,3}
>>>> In[10]:=
>>> g1[1,6]//Timing
>>>> Out[10]=
>>> {0.005429 Second,5}
>>>>> The correct answer is 3 (returned by the first method). The  
>>> faster method found two false solutions. This should not  
>>> matter if you are interested only in approximate answers (as  
>>> you seem to be) but it is worth keeping in mind.
>>>> Andrzej Kozlowski
>> I have noticed one more obvious improvement. We can replace  
>> test by:
>> test[n_] :=
>>     Mod[n, 4] =!= 3 && JacobiSymbol[n,
>>   Prime[Random[Integer, {2, 100}]]] \[CapitalEth] -1 && \
Element[Sqrt[n],
>> Integers]
>> and test1 by
>> test1[n_] := With[{w = Sqrt[N[n]]}, Mod[n, 4] =!= 3 && w ==  
>> Round[w]
>>    ]
>> We are simply making use of the easily to prove fact that an  
>> integer of the form 4 k + 1 can never be the sum of two  
>> squares. There is a noticeable improvement in the performance  
>> of g:
>> g[3,2]//Timing
>> {58.0786 Second,283}
>> However, the performance of g1 seems to actually decline  
>> slightly:
>> g1[3,2]//Timing
>> {40.8776 Second,283}
>> However, there are fewer cases of \false solutions\:
>> In[22]:=
>> g1[1,6]//Timing
>> Out[22]=
>> {0.006694 Second,4}
>> I am still not sure what is the most efficient use of  
>> JacobiSymbol in this kind of problems. In the first version of  
>> my code I used a test involving the first few odd primes,  
>> afterwards I switched to using just one random prime in taken  
>> form a certain range. Evaluation of JacobiSympol[m,p] != -1 is  
>> certainly much faster than that of Element[Sqrt[m],Integers],  
>> but it is not clear what is the most efficient number of  
>> primes to use, which are the best primes and whether it is  
>> better to choose them at random or just use a fixed selection.
>> The numerical test Sqrt[N[n]]==Round[Sqrt[N[n]] is even  
>> faster, but it will sometimes produce \false squares\.
>> Andrzej Kozlowski
>>
>
===
Subject: Re:  Re: Finding the Number of Pythagorean Triples below a bound
I need to add some corrections. There were some some small mistakes  
in the code I posted earlier, which caused problems with running the  
compiled code for negative values of k, and which also probably  
accounted for the slightly incorrect answers which the \approximate\  
code returned. I attributed the inaccuracy to the use of numerical  
precision in the test for a number being a perfect square but it  
seems that the cause was (probably) elsewhere. I don't want to try to  
make this message too long so I won't bother explaining what I think  
the mistakes were; but I will give what I think is the correct code.  
I have decided to separate the code for negative and positive values  
of k. The code for negative k works also for positive k's but is  
slightly slower, due to the extra test that needs to be performed.  
For this reason, and for the sake of greater clarity I have decided  
to separate the two codes.
The code for positive k:
countTriplesP = Compile[{{m, _Integer}, {k, _Integer}}, Module[
   {i = 0, c2, diff, sdiff},
   Do [
    If[Mod[c^2 + k, 4] != 3, c2 = Floor[Sqrt[(c^2 + k)/2]];
     Do [
       diff = c^2 + k - b^2;
          sdiff = Sqrt[N[diff]];
          If [sdiff >= b && sdiff == Round[sdiff], i++];
       , {b, 1, c2}]];
, {c, 1, 10^m - 1}];
i]]
The code for negative k:
countTriplesN = Compile[{{m, _Integer}, {k, _Integer}}, Module[
   {i = 1, c2, diff, sdiff},
   Do [
    If[Mod[c^2 + k, 4] != 3 && c^2 + k Ò 0, c2 = \
Floor[Sqrt[(c^2 + k)/ 
2]];
     Do [
       diff = c^2 + k - b^2;
          sdiff = Sqrt[N[diff]];
          If [sdiff >= b && sdiff == Round[sdiff], i++];
       , {b, 1, c2}]];
, {c, 1, 10^m - 1}];
i]]
Now we get:
countTriplesP[1,6]+countTriplesN[1,-6]//Timing
{0.000221 Second,3}
countTriplesP[3,2]+countTriplesN[3,-2]//Timing
{0.95186 Second,282}
countTriplesP[4,2]+countTriplesN[4,-2]//Timing
{95.2177 Second,2762}
Note that these values are consistently less by one form the values  
obtained by Titus and also by me using my earlier \exact\ code, but  
actually I believe that this disparity was due to mistakes in the  
earlier code. In any case, if we replace the \numerical\ test for a  
perfect square by the much slower \exact\ test, the answers will be  
the same, so the difference of 1 is certainly not due to the use of  
numerical precision. Anyway, everything works fast and looks  
perfectly satisfactory but then there is a snag. I decided to run the  
code for m=5 and k=2, started it and went out for several hours. When  
I came back I was rather disappointed to see that it was still  
running and then I saw the message:
countTriplesP[5,2]//Timing
CompiledFunction::\cfn\ Numerical error  encountered at instruction  
10; proceeding with uncompiled evaluation.
I assume the cause of this is that on 32 bit computers Compile cannot  
deal with integers larger than 2^32 but we have:
c=10^5;
Log[2,c^2]//N
33.2193
I can't at the moment see any way around this problem except to run  
uncompiled code (far too long) or try a 64 bit computer.  
Unfortunately I do not have one and I don't think Titus has one, so  
if any body has a 64-bit computer with a 64 bit version of  
Mathematica installed, and a little spare time, I think both of us  
would like to know how the above code performs in this setting. If it  
works it will provide a nice bit of advertising for 64 bit computing.
Andrzej Kozlowski
> My latest code has just managed:
 countT[4,2]+countT[4,-2]//Timing

> {93.9638 Second,2762}

> That is less than two minutes on a 1 gigahertz computer. The  
> correct answer is actually 2763 (by my earlier computation using  
> the exact test) so we have lost one solution but gained more than  
> 50 fold improvement in performance!
> The case m=5 is now certainly feasible, although I am not sure if I  
> wish my not very powerful PowerBook to be occupied for so long, as  
> I need to use Mathematica for other tasks. Perhaps I can now leave  
> this to others.
 Andrzej Kozlowski

>> it becomes really fast ;-)
>> countT = Compile[{{m, _Integer}, {k, _Integer}}, Module[
>>   {i = 0, c2, diff, sdiff},
>>   Do [
>>    If[Mod[c^2 + k, 4] != 3, c2 = Floor[Sqrt[(c^2 + k)/2]];
>>     Do [
>>       diff = c^2 + k - b^2;
>>          sdiff = Sqrt[N[diff]];
>>          If [sdiff >= b && sdiff == Round[sdiff], i++];
>>       , {b, 1, c2}]];
>> , {c, 1, 10^m - 1}];
>> i]]
>> countT[3,2]+countT[3,-2]//Timing
>> {1.0032 Second,282}
>> I will try the case m=4 and m=5 and send you the results. I am not  
>> promising to do this very soon, just in case ;-)
>> Andrzej Kozlowski
> Here is the fastest code I have seen so far that works (it is  
> based on the one that Daniel sent you). I have corrected and  
> enhanced  and enhanced it. It gives approximate answers for  
> reasons that I explained in earlier postings (use of machine  
> arithmetic to test for perfect integers). Of course it is easy to  
> replace the code by exact code by replacing the numerical test  
> for a perfect square by an exact one.

> countTriples[m_,k_] := Module[
>   {i=0, c2, diff, sdiff},
>   Do [
>    If[Mod[c^2+k,4]!=3, c2 = Floor[Sqrt[(c^2+k)/2]];
>     Do [
>       diff = c^2+k-b^2;
>          sdiff = Sqrt[N[diff]];
>          If [sdiff>=b&&sdiff==Round[sdiff],i++];
>       , {b,1,c2}]];
> ,{c,1,10^m-1}];
> i]

> countTriples[3,2]+countTriples[3,-2]//Timing

> {12.3746 Second,282}
 The correct answer is 283.
 This code should easily deal with the case m=4 (I have not yet  
> tried it) and I think even m=5 should now be within reach.
 Andrzej Kozlowski


>> The \improvement\ below which I sent a little earlier is wrong  
>> (even though it returned correct answers). Obviously the point  
>> is that c^2+k  can't be of the form 4n + 3, but there is no  
>> reason why a^2+b^2-k can't be of that form. Since my code does  
>> not explicitly select c it can't make use of this additional  
>> improvement. A different code, which uses explicit choices of  
>> (say) a and c and tests for b being a perfect square could  
>> exploit this fact and perhaps gain extra speed. It should not be  
>> difficult to write such a code along the lines I have been using.
>> Andrzej
>>>>>>>>>> Hello all,
>>>> emailed
>>> me.
>>>> To recall, the problem was counting the number of solutions  
>>> to a
>>> bivariate polynomial equal to a square,
>>>> Poly(a,b) = c^2
>>>> One form that interested me was the Pythagorean-like equation:
>>>> a^2 + b^2 = c^2 + k
>>>> for {a,b} a positive integer, 0<a<=b, and k any small  
>>> integer. I was
>>> wondering about the density of solutions to this since I  
>>> knew in the
>>> special case of k=0, let S(N) be the number of primitive  
>>> solutions
>>> with
>>> c < N, then S(N)/N = 1/(2pi) as N -> inf.
>>>> For k a squarefree integer, it is convenient that any  
>>> solution is
>>> also
>>> primitive. I used a simple code that allowed me to find S 
>>> (10^m) with
>>> m=1,2,3 for small values of k (for m=4 took my code more  
>>> than 30 mins
>>> so I aborted it). The data is given below:
>>>> Note: Values are total S(N) for *both* k & -k:
>>>> k = 2
>>> S(N) = 4, 30, 283
>>>> k = 3
>>> S(N) = 3, 41, 410
>>>> k = 5
>>> S(N) = 3, 43, 426
>>>> k = 6
>>> S(N) = 3, 36, 351
>>>> Question: Does S(N)/N for these also converge? For example,  
>>> for the
>>> particular case of k = -6, we have
>>>> S(N) = 2, 20, 202
>>>> which looks suspiciously like the ratio might be converging.
>>>> Anybody know of a code for this that can find m=4,5,6 in a  
>>> reasonable
>>> amount of time?
>>>>> Yours,
>>>> Titus
>>>>>>>>> Here is a piece code which utilises the ideas I have  
>>>> described in
>>>> my previous posts:
>>>>>> ls = Prime /@ Range[3, 10];
>>>>>> test[n_] :=
>>>>       Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>>> Integers]
>>>>>> f[P_, k_] := Sum[If[(w =
>>>>    a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a,
>>>>           Floor[Sqrt[P^2 - a^2]]}]
>>>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>>>> We can easily confirm the results of your computations,.e.g.  
>>>> for k=2.
>>>>>>>> Table[g[i,2],{i,3}]
>>>>>>>> {4,30,283}
>>>>>>>> Since you have not revealed the \simple code\ you have used  
>>>> it is
>>>> hard to tell if the above one is any better. It is however,
>>>> certainly capable of solving the problem for m=4:
>>>>>>>> g[4,2]//Timing
>>>>>>>> {4779.39 Second,2763}
>>>>>> The time it took on my 1 Gigahertz PowerBook was over 70  
>>>> minutes,
>>>> which is longer than you thought \reasonable\, so I am still  
>>>> not
>>>> sure if this is any improvement on what you already have.  
>>>> The time
>>>> complexity of this algorithm seems somewhat larger than  
>>>> exponential
>>>> so I would expect that it will take about 6 hours to deal  
>>>> with n=5
>>>> on my computer, and perhaps 2 weeks to deal with n=6.
>>>>>> Andrzej Kozlowski
>>>>>>>>>>> I mistakenly copied and pasted a wrong (earlier) definition  
>>> of f.
>>> Here is the correct one:
>>>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>      1 && (w = a^2 +
>>>           b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>>            P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]]
>>>> The definition of g is as before:
>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>> Andrzej Kozlowski
>>> Below is a faster version of the above code. (It owes a  
>> significant improvement to Daniel Lichtblau, which I borrowed  
>> from him without his knowledge ;-))
>> test[n_] :=
>>     JacobiSymbol[n,
>>   Prime[Random[Integer, {2, 20}]]] \[CapitalEth] -1 && \
Element[Sqrt[n],
>> Integers]
>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>     1 && (w = a^2 +
>>          b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>> +k]]}]]
>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>> The improvement is in the upper bound on a in the sum. Since a  
>> is the smaller of the two squares whose sum is equal to P^2+k  
>> it can't be larger than  Floor[Sqrt[(P^2+k)/2]].
>> Note that you can improve the performance by loosing some  
>> accuracy if you use a cruder test for a perfect square:
>> test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w]
>>    ]
>> f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>     1 && (w = a^2 +
>>          b^2 - k) < P^2 && test1[w], 1, 0], {a, 1,
>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>> +k]]}]]
>> g1[m_, k_] := f1[10^m, k] + f1[10^m, -k]
>> Let's compare the two cases.
>> In[7]:=
>> g[3,2]//Timing
>> Out[7]=
>> {89.554 Second,283}
>> In[8]:=
>> g1[3,2]//Timing
>> Out[8]=
>> {37.376 Second,283}
>> So we see that we get the same answer and the second approach  
>> is considerably faster. However:
>> In[9]:=
>> g[1,6]//Timing
>> Out[9]=
>> {0.008863 Second,3}
>> In[10]:=
>> g1[1,6]//Timing
>> Out[10]=
>> {0.005429 Second,5}
>> The correct answer is 3 (returned by the first method). The  
>> faster method found two false solutions. This should not  
>> matter if you are interested only in approximate answers (as  
>> you seem to be) but it is worth keeping in mind.
>> Andrzej Kozlowski
>>>>> I have noticed one more obvious improvement. We can replace  
>>> test by:
>>>> test[n_] :=
>>>     Mod[n, 4] =!= 3 && JacobiSymbol[n,
>>>   Prime[Random[Integer, {2, 100}]]] \[CapitalEth] -1 && \
Element[Sqrt[n],
>>> Integers]
>>>> and test1 by
>>>> test1[n_] := With[{w = Sqrt[N[n]]}, Mod[n, 4] =!= 3 && w ==  
>>> Round[w]
>>>    ]
>>>> We are simply making use of the easily to prove fact that an  
>>> integer of the form 4 k + 1 can never be the sum of two  
>>> squares. There is a noticeable improvement in the performance  
>>> of g:
>>>>> g[3,2]//Timing
>>>> {58.0786 Second,283}
>>>> However, the performance of g1 seems to actually decline slightly:
>>>>> g1[3,2]//Timing
>>>> {40.8776 Second,283}
>>>>> However, there are fewer cases of \false solutions\:
>>>> In[22]:=
>>> g1[1,6]//Timing
>>>> Out[22]=
>>> {0.006694 Second,4}
>>>>> I am still not sure what is the most efficient use of  
>>> JacobiSymbol in this kind of problems. In the first version of  
>>> my code I used a test involving the first few odd primes,  
>>> afterwards I switched to using just one random prime in taken  
>>> form a certain range. Evaluation of JacobiSympol[m,p] != -1 is  
>>> certainly much faster than that of Element[Sqrt[m],Integers],  
>>> but it is not clear what is the most efficient number of primes  
>>> to use, which are the best primes and whether it is better to  
>>> choose them at random or just use a fixed selection.
>>> The numerical test Sqrt[N[n]]==Round[Sqrt[N[n]] is even faster,  
>>> but it will sometimes produce \false squares\.
>>>> Andrzej Kozlowski
>>>
>
===
Subject: Re:  Re: Finding the Number of Pythagorean Triples below a bound
My latest code has just managed:
countT[4,2]+countT[4,-2]//Timing
{93.9638 Second,2762}
That is less than two minutes on a 1 gigahertz computer. The correct  
answer is actually 2763 (by my earlier computation using the exact  
test) so we have lost one solution but gained more than 50 fold  
improvement in performance!
The case m=5 is now certainly feasible, although I am not sure if I  
wish my not very powerful PowerBook to be occupied for so long, as I  
need to use Mathematica for other tasks. Perhaps I can now leave this  
to others.
Andrzej Kozlowski
> it becomes really fast ;-)
 countT = Compile[{{m, _Integer}, {k, _Integer}}, Module[
>   {i = 0, c2, diff, sdiff},
>   Do [
>    If[Mod[c^2 + k, 4] != 3, c2 = Floor[Sqrt[(c^2 + k)/2]];
>     Do [
>       diff = c^2 + k - b^2;
>          sdiff = Sqrt[N[diff]];
>          If [sdiff >= b && sdiff == Round[sdiff], i++];
>       , {b, 1, c2}]];
> , {c, 1, 10^m - 1}];
> i]]

> countT[3,2]+countT[3,-2]//Timing

> {1.0032 Second,282}
 I will try the case m=4 and m=5 and send you the results. I am not  
> promising to do this very soon, just in case ;-)
 Andrzej Kozlowski


>> Here is the fastest code I have seen so far that works (it is  
>> based on the one that Daniel sent you). I have corrected and  
>> enhanced  and enhanced it. It gives approximate answers for  
>> reasons that I explained in earlier postings (use of machine  
>> arithmetic to test for perfect integers). Of course it is easy to  
>> replace the code by exact code by replacing the numerical test for  
>> a perfect square by an exact one.
>> countTriples[m_,k_] := Module[
>>   {i=0, c2, diff, sdiff},
>>   Do [
>>    If[Mod[c^2+k,4]!=3, c2 = Floor[Sqrt[(c^2+k)/2]];
>>     Do [
>>       diff = c^2+k-b^2;
>>          sdiff = Sqrt[N[diff]];
>>          If [sdiff>=b&&sdiff==Round[sdiff],i++];
>>       , {b,1,c2}]];
>> ,{c,1,10^m-1}];
>> i]
>> countTriples[3,2]+countTriples[3,-2]//Timing
>> {12.3746 Second,282}
>> The correct answer is 283.
>> This code should easily deal with the case m=4 (I have not yet  
>> tried it) and I think even m=5 should now be within reach.
>> Andrzej Kozlowski
> The \improvement\ below which I sent a little earlier is wrong  
> (even though it returned correct answers). Obviously the point is  
> that c^2+k  can't be of the form 4n + 3, but there is no reason  
> why a^2+b^2-k can't be of that form. Since my code does not  
> explicitly select c it can't make use of this additional  
> improvement. A different code, which uses explicit choices of  
> (say) a and c and tests for b being a perfect square could  
> exploit this fact and perhaps gain extra speed. It should not be  
> difficult to write such a code along the lines I have been using.
 Andrzej

>>>>>>> Hello all,
>>>>>> emailed
>>>> me.
>>>>>> To recall, the problem was counting the number of solutions  
>>>> to a
>>>> bivariate polynomial equal to a square,
>>>>>> Poly(a,b) = c^2
>>>>>> One form that interested me was the Pythagorean-like equation:
>>>>>> a^2 + b^2 = c^2 + k
>>>>>> for {a,b} a positive integer, 0<a<=b, and k any small  
>>>> integer. I was
>>>> wondering about the density of solutions to this since I  
>>>> knew in the
>>>> special case of k=0, let S(N) be the number of primitive  
>>>> solutions
>>>> with
>>>> c < N, then S(N)/N = 1/(2pi) as N -> inf.
>>>>>> For k a squarefree integer, it is convenient that any  
>>>> solution is
>>>> also
>>>> primitive. I used a simple code that allowed me to find S 
>>>> (10^m) with
>>>> m=1,2,3 for small values of k (for m=4 took my code more  
>>>> than 30 mins
>>>> so I aborted it). The data is given below:
>>>>>> Note: Values are total S(N) for *both* k & -k:
>>>>>> k = 2
>>>> S(N) = 4, 30, 283
>>>>>> k = 3
>>>> S(N) = 3, 41, 410
>>>>>> k = 5
>>>> S(N) = 3, 43, 426
>>>>>> k = 6
>>>> S(N) = 3, 36, 351
>>>>>> Question: Does S(N)/N for these also converge? For example,  
>>>> for the
>>>> particular case of k = -6, we have
>>>>>> S(N) = 2, 20, 202
>>>>>> which looks suspiciously like the ratio might be converging.
>>>>>> Anybody know of a code for this that can find m=4,5,6 in a  
>>>> reasonable
>>>> amount of time?
>>>>>>>> Yours,
>>>>>> Titus
>>>>>>> Here is a piece code which utilises the ideas I have  
>>> described in
>>> my previous posts:
>>>> ls = Prime /@ Range[3, 10];
>>>> test[n_] :=
>>>       Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>> Integers]
>>>> f[P_, k_] := Sum[If[(w =
>>>    a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a,
>>>           Floor[Sqrt[P^2 - a^2]]}]
>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>> We can easily confirm the results of your computations,.e.g.  
>>> for k=2.
>>>>> Table[g[i,2],{i,3}]
>>>>> {4,30,283}
>>>>> Since you have not revealed the \simple code\ you have used  
>>> it is
>>> hard to tell if the above one is any better. It is however,
>>> certainly capable of solving the problem for m=4:
>>>>> g[4,2]//Timing
>>>>> {4779.39 Second,2763}
>>>> The time it took on my 1 Gigahertz PowerBook was over 70  
>>> minutes,
>>> which is longer than you thought \reasonable\, so I am still not
>>> sure if this is any improvement on what you already have. The  
>>> time
>>> complexity of this algorithm seems somewhat larger than  
>>> exponential
>>> so I would expect that it will take about 6 hours to deal  
>>> with n=5
>>> on my computer, and perhaps 2 weeks to deal with n=6.
>>>> Andrzej Kozlowski
>>>>> I mistakenly copied and pasted a wrong (earlier) definition of f.
>> Here is the correct one:
>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>      1 && (w = a^2 +
>>           b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>            P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]]
>> The definition of g is as before:
>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>> Andrzej Kozlowski
>>>>> Below is a faster version of the above code. (It owes a  
>>> significant improvement to Daniel Lichtblau, which I borrowed  
>>> from him without his knowledge ;-))
>>>> test[n_] :=
>>>     JacobiSymbol[n,
>>>   Prime[Random[Integer, {2, 20}]]] \[CapitalEth] -1 && \
Element[Sqrt[n],
>>> Integers]
>>>>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>     1 && (w = a^2 +
>>>          b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>>> +k]]}]]
>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>> The improvement is in the upper bound on a in the sum. Since a  
>>> is the smaller of the two squares whose sum is equal to P^2+k  
>>> it can't be larger than  Floor[Sqrt[(P^2+k)/2]].
>>>> Note that you can improve the performance by loosing some  
>>> accuracy if you use a cruder test for a perfect square:
>>>> test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w]
>>>    ]
>>>> f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>>     1 && (w = a^2 +
>>>          b^2 - k) < P^2 && test1[w], 1, 0], {a, 1,
>>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
>>> +k]]}]]
>>>>> g1[m_, k_] := f1[10^m, k] + f1[10^m, -k]
>>>>> Let's compare the two cases.
>>>> In[7]:=
>>> g[3,2]//Timing
>>>> Out[7]=
>>> {89.554 Second,283}
>>>> In[8]:=
>>> g1[3,2]//Timing
>>>> Out[8]=
>>> {37.376 Second,283}
>>>> So we see that we get the same answer and the second approach  
>>> is considerably faster. However:
>>>> In[9]:=
>>> g[1,6]//Timing
>>>> Out[9]=
>>> {0.008863 Second,3}
>>>> In[10]:=
>>> g1[1,6]//Timing
>>>> Out[10]=
>>> {0.005429 Second,5}
>>>>> The correct answer is 3 (returned by the first method). The  
>>> faster method found two false solutions. This should not matter  
>>> if you are interested only in approximate answers (as you seem  
>>> to be) but it is worth keeping in mind.
>>>> Andrzej Kozlowski
>> I have noticed one more obvious improvement. We can replace test  
>> by:
>> test[n_] :=
>>     Mod[n, 4] =!= 3 && JacobiSymbol[n,
>>   Prime[Random[Integer, {2, 100}]]] \[CapitalEth] -1 && \
Element[Sqrt[n],
>> Integers]
>> and test1 by
>> test1[n_] := With[{w = Sqrt[N[n]]}, Mod[n, 4] =!= 3 && w == Round 
>> [w]
>>    ]
>> We are simply making use of the easily to prove fact that an  
>> integer of the form 4 k + 1 can never be the sum of two squares.  
>> There is a noticeable improvement in the performance of g:
>> g[3,2]//Timing
>> {58.0786 Second,283}
>> However, the performance of g1 seems to actually decline slightly:
>> g1[3,2]//Timing
>> {40.8776 Second,283}
>> However, there are fewer cases of \false solutions\:
>> In[22]:=
>> g1[1,6]//Timing
>> Out[22]=
>> {0.006694 Second,4}
>> I am still not sure what is the most efficient use of  
>> JacobiSymbol in this kind of problems. In the first version of  
>> my code I used a test involving the first few odd primes,  
>> afterwards I switched to using just one random prime in taken  
>> form a certain range. Evaluation of JacobiSympol[m,p] != -1 is  
>> certainly much faster than that of Element[Sqrt[m],Integers],  
>> but it is not clear what is the most efficient number of primes  
>> to use, which are the best primes and whether it is better to  
>> choose them at random or just use a fixed selection.
>> The numerical test Sqrt[N[n]]==Round[Sqrt[N[n]] is even faster,  
>> but it will sometimes produce \false squares\.
>> Andrzej Kozlowski
>
>
===
Subject: Re:  Re: Finding the Number of Pythagorean Triples below a bound
Here is the fastest code I have seen so far that works (it is based  
on the one that Daniel sent you). I have corrected and enhanced  and  
enhanced it. It gives approximate answers for reasons that I  
explained in earlier postings (use of machine arithmetic to test for  
perfect integers). Of course it is easy to replace the code by exact  
code by replacing the numerical test for a perfect square by an exact  
one.
countTriples[m_,k_] := Module[
   {i=0, c2, diff, sdiff},
   Do [
    If[Mod[c^2+k,4]!=3, c2 = Floor[Sqrt[(c^2+k)/2]];
     Do [
       diff = c^2+k-b^2;
          sdiff = Sqrt[N[diff]];
          If [sdiff>=b&&sdiff==Round[sdiff],i++];
       , {b,1,c2}]];
,{c,1,10^m-1}];
i]
countTriples[3,2]+countTriples[3,-2]//Timing
{12.3746 Second,282}
The correct answer is 283.
This code should easily deal with the case m=4 (I have not yet tried  
it) and I think even m=5 should now be within reach.
Andrzej Kozlowski
> The \improvement\ below which I sent a little earlier is wrong  
> (even though it returned correct answers). Obviously the point is  
> that c^2+k  can't be of the form 4n + 3, but there is no reason why  
> a^2+b^2-k can't be of that form. Since my code does not explicitly  
> select c it can't make use of this additional improvement. A  
> different code, which uses explicit choices of (say) a and c and  
> tests for b being a perfect square could exploit this fact and  
> perhaps gain extra speed. It should not be difficult to write such  
> a code along the lines I have been using.
 Andrzej


>>> Hello all,
>> emailed
>> me.
>> To recall, the problem was counting the number of solutions to a
>> bivariate polynomial equal to a square,
>> Poly(a,b) = c^2
>> One form that interested me was the Pythagorean-like equation:
>> a^2 + b^2 = c^2 + k
>> for {a,b} a positive integer, 0<a<=b, and k any small integer.  
>> I was
>> wondering about the density of solutions to this since I knew  
>> in the
>> special case of k=0, let S(N) be the number of primitive  
>> solutions
>> with
>> c < N, then S(N)/N = 1/(2pi) as N -> inf.
>> For k a squarefree integer, it is convenient that any solution is
>> also
>> primitive. I used a simple code that allowed me to find S 
>> (10^m) with
>> m=1,2,3 for small values of k (for m=4 took my code more than  
>> 30 mins
>> so I aborted it). The data is given below:
>> Note: Values are total S(N) for *both* k & -k:
>> k = 2
>> S(N) = 4, 30, 283
>> k = 3
>> S(N) = 3, 41, 410
>> k = 5
>> S(N) = 3, 43, 426
>> k = 6
>> S(N) = 3, 36, 351
>> Question: Does S(N)/N for these also converge? For example,  
>> for the
>> particular case of k = -6, we have
>> S(N) = 2, 20, 202
>> which looks suspiciously like the ratio might be converging.
>> Anybody know of a code for this that can find m=4,5,6 in a  
>> reasonable
>> amount of time?
>> Yours,
>> Titus
>>>>> Here is a piece code which utilises the ideas I have described in
>>> my previous posts:
>>>> ls = Prime /@ Range[3, 10];
>>>> test[n_] :=
>>>       Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>>> Integers]
>>>> f[P_, k_] := Sum[If[(w =
>>>    a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a,
>>>           Floor[Sqrt[P^2 - a^2]]}]
>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>>>> We can easily confirm the results of your computations,.e.g.  
>>> for k=2.
>>>>> Table[g[i,2],{i,3}]
>>>>> {4,30,283}
>>>>> Since you have not revealed the \simple code\ you have used it is
>>> hard to tell if the above one is any better. It is however,
>>> certainly capable of solving the problem for m=4:
>>>>> g[4,2]//Timing
>>>>> {4779.39 Second,2763}
>>>> The time it took on my 1 Gigahertz PowerBook was over 70 minutes,
>>> which is longer than you thought \reasonable\, so I am still not
>>> sure if this is any improvement on what you already have. The time
>>> complexity of this algorithm seems somewhat larger than  
>>> exponential
>>> so I would expect that it will take about 6 hours to deal with n=5
>>> on my computer, and perhaps 2 weeks to deal with n=6.
>>>> Andrzej Kozlowski
>>>>> I mistakenly copied and pasted a wrong (earlier) definition of f.
>> Here is the correct one:
>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>      1 && (w = a^2 +
>>           b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>            P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]]
>> The definition of g is as before:
>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>> Andrzej Kozlowski
>
 Below is a faster version of the above code. (It owes a  
> significant improvement to Daniel Lichtblau, which I borrowed  
> from him without his knowledge ;-))
 test[n_] :=
>     JacobiSymbol[n,
>   Prime[Random[Integer, {2, 20}]]] \[CapitalEth] -1 && Element[Sqrt[n],
> Integers]

> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>     1 && (w = a^2 +
>          b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
> +k]]}]]
 g[m_, k_] := f[10^m, k] + f[10^m, -k]
 The improvement is in the upper bound on a in the sum. Since a is  
> the smaller of the two squares whose sum is equal to P^2+k it  
> can't be larger than  Floor[Sqrt[(P^2+k)/2]].
 Note that you can improve the performance by loosing some  
> accuracy if you use a cruder test for a perfect square:
 test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w]
>    ]
 f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>     1 && (w = a^2 +
>          b^2 - k) < P^2 && test1[w], 1, 0], {a, 1,
>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 
> +k]]}]]

> g1[m_, k_] := f1[10^m, k] + f1[10^m, -k]

> Let's compare the two cases.
 In[7]:=
> g[3,2]//Timing
 Out[7]=
> {89.554 Second,283}
 In[8]:=
> g1[3,2]//Timing
 Out[8]=
> {37.376 Second,283}
 So we see that we get the same answer and the second approach is  
> considerably faster. However:
 In[9]:=
> g[1,6]//Timing
 Out[9]=
> {0.008863 Second,3}
 In[10]:=
> g1[1,6]//Timing
 Out[10]=
> {0.005429 Second,5}

> The correct answer is 3 (returned by the first method). The  
> faster method found two false solutions. This should not matter  
> if you are interested only in approximate answers (as you seem to  
> be) but it is worth keeping in mind.
 Andrzej Kozlowski
>> I have noticed one more obvious improvement. We can replace test by:
>> test[n_] :=
>>     Mod[n, 4] =!= 3 && JacobiSymbol[n,
>>   Prime[Random[Integer, {2, 100}]]] \[CapitalEth] -1 && Element[Sqrt[n],
>> Integers]
>> and test1 by
>> test1[n_] := With[{w = Sqrt[N[n]]}, Mod[n, 4] =!= 3 && w == Round[w]
>>    ]
>> We are simply making use of the easily to prove fact that an  
>> integer of the form 4 k + 1 can never be the sum of two squares.  
>> There is a noticeable improvement in the performance of g:
>> g[3,2]//Timing
>> {58.0786 Second,283}
>> However, the performance of g1 seems to actually decline slightly:
>> g1[3,2]//Timing
>> {40.8776 Second,283}
>> However, there are fewer cases of \false solutions\:
>> In[22]:=
>> g1[1,6]//Timing
>> Out[22]=
>> {0.006694 Second,4}
>> I am still not sure what is the most efficient use of JacobiSymbol  
>> in this kind of problems. In the first version of my code I used a  
>> test involving the first few odd primes, afterwards I switched to  
>> using just one random prime in taken form a certain range.  
>> Evaluation of JacobiSympol[m,p] != -1 is certainly much faster  
>> than that of Element[Sqrt[m],Integers], but it is not clear what  
>> is the most efficient number of primes to use, which are the best  
>> primes and whether it is better to choose them at random or just  
>> use a fixed selection.
>> The numerical test Sqrt[N[n]]==Round[Sqrt[N[n]] is even faster,  
>> but it will sometimes produce \false squares\.
>> Andrzej Kozlowski
>
===
Subject: Re:  Re: Finding the Number of Pythagorean Triples below a bound
> Hello Andrzej,
 The code I used was this:
 c1=Sqrt[N[a^2+b^2-k]];
 countTriples[m_]:=Block[{a,b,i=0},For[a=1,a<(10^m/Sqrt[2]),a++,
> For[b=a,b<10^m,b++,If[c1==Round[c1]&&c1<10^m,i++]]];i]
 though m=3 takes about 50 sec (in an old comp) while m=4 will  
> prolly be
 around 5000 sec (83 mins).

> -Titus
>
This is very odd indeed. For example, I get:
c1=Sqrt[N[a^2+b^2-k]];
countTriples[m_]:=Block[{a,b,i=0},For[a=1,a<(10^m/Sqrt[2]),a++,
For[b=a,b<10^m,b++,If[c1==Round[c1]&&c1<10^m,i++]]];i]
In[6]:=
k=3;
countTriples[2]
Invalid comparison
     with  0. I attempted.
13
The correct answer is:
f[10^2,-3]
28
Similarly for other values. Your code usually doesn't work on my  
machine and when it does it produces wrong answers.
Andrzej Kozlowski
===
Subject: Re:  Re: Finding the Number of Pythagorean Triples below a bound
The \improvement\ below which I sent a little earlier is wrong (even  
though it returned correct answers). Obviously the point is that c^2 
+k  can't be of the form 4n + 3, but there is no reason why a^2+b^2-k  
can't be of that form. Since my code does not explicitly select c it  
can't make use of this additional improvement. A different code,  
which uses explicit choices of (say) a and c and tests for b being a  
perfect square could exploit this fact and perhaps gain extra speed.  
It should not be difficult to write such a code along the lines I  
have been using.
Andrzej

>
>>> Hello all,
>>>> emailed
>>> me.
>>>> To recall, the problem was counting the number of solutions to a
>>> bivariate polynomial equal to a square,
>>>> Poly(a,b) = c^2
>>>> One form that interested me was the Pythagorean-like equation:
>>>> a^2 + b^2 = c^2 + k
>>>> for {a,b} a positive integer, 0<a<=b, and k any small integer.  
>>> I was
>>> wondering about the density of solutions to this since I knew  
>>> in the
>>> special case of k=0, let S(N) be the number of primitive solutions
>>> with
>>> c < N, then S(N)/N = 1/(2pi) as N -> inf.
>>>> For k a squarefree integer, it is convenient that any solution is
>>> also
>>> primitive. I used a simple code that allowed me to find S(10^m)  
>>> with
>>> m=1,2,3 for small values of k (for m=4 took my code more than  
>>> 30 mins
>>> so I aborted it). The data is given below:
>>>> Note: Values are total S(N) for *both* k & -k:
>>>> k = 2
>>> S(N) = 4, 30, 283
>>>> k = 3
>>> S(N) = 3, 41, 410
>>>> k = 5
>>> S(N) = 3, 43, 426
>>>> k = 6
>>> S(N) = 3, 36, 351
>>>> Question: Does S(N)/N for these also converge? For example, for  
>>> the
>>> particular case of k = -6, we have
>>>> S(N) = 2, 20, 202
>>>> which looks suspiciously like the ratio might be converging.
>>>> Anybody know of a code for this that can find m=4,5,6 in a  
>>> reasonable
>>> amount of time?
>>>>> Yours,
>>>> Titus
>>> Here is a piece code which utilises the ideas I have described in
>> my previous posts:
>> ls = Prime /@ Range[3, 10];
>> test[n_] :=
>>       Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>> Integers]
>> f[P_, k_] := Sum[If[(w =
>>    a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a,
>>           Floor[Sqrt[P^2 - a^2]]}]
>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>> We can easily confirm the results of your computations,.e.g. for  
>> k=2.
>> Table[g[i,2],{i,3}]
>> {4,30,283}
>> Since you have not revealed the \simple code\ you have used it is
>> hard to tell if the above one is any better. It is however,
>> certainly capable of solving the problem for m=4:
>> g[4,2]//Timing
>> {4779.39 Second,2763}
>> The time it took on my 1 Gigahertz PowerBook was over 70 minutes,
>> which is longer than you thought \reasonable\, so I am still not
>> sure if this is any improvement on what you already have. The time
>> complexity of this algorithm seems somewhat larger than exponential
>> so I would expect that it will take about 6 hours to deal with n=5
>> on my computer, and perhaps 2 weeks to deal with n=6.
>> Andrzej Kozlowski
>
 I mistakenly copied and pasted a wrong (earlier) definition of f.
> Here is the correct one:
 f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>      1 && (w = a^2 +
>           b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>            P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]]
 The definition of g is as before:
 g[m_, k_] := f[10^m, k] + f[10^m, -k]
 Andrzej Kozlowski
> Below is a faster version of the above code. (It owes a  
>> significant improvement to Daniel Lichtblau, which I borrowed from  
>> him without his knowledge ;-))
>> test[n_] :=
>>     JacobiSymbol[n,
>>   Prime[Random[Integer, {2, 20}]]] \[CapitalEth] -1 && Element[Sqrt[n],
>> Integers]
>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>     1 && (w = a^2 +
>>          b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2+k]]}]]
>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>> The improvement is in the upper bound on a in the sum. Since a is  
>> the smaller of the two squares whose sum is equal to P^2+k it  
>> can't be larger than  Floor[Sqrt[(P^2+k)/2]].
>> Note that you can improve the performance by loosing some accuracy  
>> if you use a cruder test for a perfect square:
>> test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w]
>>    ]
>> f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>     1 && (w = a^2 +
>>          b^2 - k) < P^2 && test1[w], 1, 0], {a, 1,
>>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2+k]]}]]
>> g1[m_, k_] := f1[10^m, k] + f1[10^m, -k]
>> Let's compare the two cases.
>> In[7]:=
>> g[3,2]//Timing
>> Out[7]=
>> {89.554 Second,283}
>> In[8]:=
>> g1[3,2]//Timing
>> Out[8]=
>> {37.376 Second,283}
>> So we see that we get the same answer and the second approach is  
>> considerably faster. However:
>> In[9]:=
>> g[1,6]//Timing
>> Out[9]=
>> {0.008863 Second,3}
>> In[10]:=
>> g1[1,6]//Timing
>> Out[10]=
>> {0.005429 Second,5}
>> The correct answer is 3 (returned by the first method). The faster  
>> method found two false solutions. This should not matter if you  
>> are interested only in approximate answers (as you seem to be) but  
>> it is worth keeping in mind.
>> Andrzej Kozlowski

> I have noticed one more obvious improvement. We can replace test by:
 test[n_] :=
>     Mod[n, 4] =!= 3 && JacobiSymbol[n,
>   Prime[Random[Integer, {2, 100}]]] \[CapitalEth] -1 && Element[Sqrt[n],
> Integers]
 and test1 by
 test1[n_] := With[{w = Sqrt[N[n]]}, Mod[n, 4] =!= 3 && w == Round[w]
>    ]
 We are simply making use of the easily to prove fact that an  
> integer of the form 4 k + 1 can never be the sum of two squares.  
> There is a noticeable improvement in the performance of g:

> g[3,2]//Timing
 {58.0786 Second,283}
 However, the performance of g1 seems to actually decline slightly:

> g1[3,2]//Timing
 {40.8776 Second,283}

> However, there are fewer cases of \false solutions\:
 In[22]:=
> g1[1,6]//Timing
 Out[22]=
> {0.006694 Second,4}

> I am still not sure what is the most efficient use of JacobiSymbol  
> in this kind of problems. In the first version of my code I used a  
> test involving the first few odd primes, afterwards I switched to  
> using just one random prime in taken form a certain range.  
> Evaluation of JacobiSympol[m,p] != -1 is certainly much faster than  
> that of Element[Sqrt[m],Integers], but it is not clear what is the  
> most efficient number of primes to use, which are the best primes  
> and whether it is better to choose them at random or just use a  
> fixed selection.
> The numerical test Sqrt[N[n]]==Round[Sqrt[N[n]] is even faster, but  
> it will sometimes produce \false squares\.
 Andrzej Kozlowski

===
Subject: Re:  Re: Finding the Number of Pythagorean Triples below a bound

>
> Hello all,
>> emailed
>> me.
>> To recall, the problem was counting the number of solutions to a
>> bivariate polynomial equal to a square,
>> Poly(a,b) = c^2
>> One form that interested me was the Pythagorean-like equation:
>> a^2 + b^2 = c^2 + k
>> for {a,b} a positive integer, 0<a<=b, and k any small integer. I  
>> was
>> wondering about the density of solutions to this since I knew in  
>> the
>> special case of k=0, let S(N) be the number of primitive solutions
>> with
>> c < N, then S(N)/N = 1/(2pi) as N -> inf.
>> For k a squarefree integer, it is convenient that any solution is
>> also
>> primitive. I used a simple code that allowed me to find S(10^m)  
>> with
>> m=1,2,3 for small values of k (for m=4 took my code more than 30  
>> mins
>> so I aborted it). The data is given below:
>> Note: Values are total S(N) for *both* k & -k:
>> k = 2
>> S(N) = 4, 30, 283
>> k = 3
>> S(N) = 3, 41, 410
>> k = 5
>> S(N) = 3, 43, 426
>> k = 6
>> S(N) = 3, 36, 351
>> Question: Does S(N)/N for these also converge? For example, for the
>> particular case of k = -6, we have
>> S(N) = 2, 20, 202
>> which looks suspiciously like the ratio might be converging.
>> Anybody know of a code for this that can find m=4,5,6 in a  
>> reasonable
>> amount of time?
>> Yours,
>> Titus
>
 Here is a piece code which utilises the ideas I have described in
> my previous posts:
 ls = Prime /@ Range[3, 10];
 test[n_] :=
>       Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
> Integers]
 f[P_, k_] := Sum[If[(w =
>    a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a,
>           Floor[Sqrt[P^2 - a^2]]}]
 g[m_, k_] := f[10^m, k] + f[10^m, -k]
 We can easily confirm the results of your computations,.e.g. for  
> k=2.

> Table[g[i,2],{i,3}]

> {4,30,283}

> Since you have not revealed the \simple code\ you have used it is
> hard to tell if the above one is any better. It is however,
> certainly capable of solving the problem for m=4:

> g[4,2]//Timing

> {4779.39 Second,2763}
 The time it took on my 1 Gigahertz PowerBook was over 70 minutes,
> which is longer than you thought \reasonable\, so I am still not
> sure if this is any improvement on what you already have. The time
> complexity of this algorithm seems somewhat larger than exponential
> so I would expect that it will take about 6 hours to deal with n=5
> on my computer, and perhaps 2 weeks to deal with n=6.
 Andrzej Kozlowski

> I mistakenly copied and pasted a wrong (earlier) definition of f.
>> Here is the correct one:
>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>>      1 && (w = a^2 +
>>           b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>>            P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]]
>> The definition of g is as before:
>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>> Andrzej Kozlowski
>
 Below is a faster version of the above code. (It owes a significant  
> improvement to Daniel Lichtblau, which I borrowed from him without  
> his knowledge ;-))
 test[n_] :=
>     JacobiSymbol[n,
>   Prime[Random[Integer, {2, 20}]]] \[CapitalEth] -1 && Element[Sqrt[n],
> Integers]

> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>     1 && (w = a^2 +
>          b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2+k]]}]]
 g[m_, k_] := f[10^m, k] + f[10^m, -k]
 The improvement is in the upper bound on a in the sum. Since a is  
> the smaller of the two squares whose sum is equal to P^2+k it can't  
> be larger than  Floor[Sqrt[(P^2+k)/2]].
 Note that you can improve the performance by loosing some accuracy  
> if you use a cruder test for a perfect square:
 test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w]
>    ]
 f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>     1 && (w = a^2 +
>          b^2 - k) < P^2 && test1[w], 1, 0], {a, 1,
>           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2+k]]}]]

> g1[m_, k_] := f1[10^m, k] + f1[10^m, -k]

> Let's compare the two cases.
 In[7]:=
> g[3,2]//Timing
 Out[7]=
> {89.554 Second,283}
 In[8]:=
> g1[3,2]//Timing
 Out[8]=
> {37.376 Second,283}
 So we see that we get the same answer and the second approach is  
> considerably faster. However:
 In[9]:=
> g[1,6]//Timing
 Out[9]=
> {0.008863 Second,3}
 In[10]:=
> g1[1,6]//Timing
 Out[10]=
> {0.005429 Second,5}

> The correct answer is 3 (returned by the first method). The faster  
> method found two false solutions. This should not matter if you are  
> interested only in approximate answers (as you seem to be) but it  
> is worth keeping in mind.
 Andrzej Kozlowski
I have noticed one more obvious improvement. We can replace test by:
test[n_] :=
     Mod[n, 4] =!= 3 && JacobiSymbol[n,
   Prime[Random[Integer, {2, 100}]]] \[CapitalEth] -1 && Element[Sqrt[n],
Integers]
and test1 by
test1[n_] := With[{w = Sqrt[N[n]]}, Mod[n, 4] =!= 3 && w == Round[w]
    ]
We are simply making use of the easily to prove fact that an integer  
of the form 4 k + 1 can never be the sum of two squares. There is a  
noticeable improvement in the performance of g:
g[3,2]//Timing
{58.0786 Second,283}
However, the performance of g1 seems to actually decline slightly:
g1[3,2]//Timing
{40.8776 Second,283}
However, there are fewer cases of \false solutions\:
In[22]:=
g1[1,6]//Timing
Out[22]=
{0.006694 Second,4}
I am still not sure what is the most efficient use of JacobiSymbol in  
this kind of problems. In the first version of my code I used a test  
involving the first few odd primes, afterwards I switched to using  
just one random prime in taken form a certain range. Evaluation of  
JacobiSympol[m,p] != -1 is certainly much faster than that of Element 
[Sqrt[m],Integers], but it is not clear what is the most efficient  
number of primes to use, which are the best primes and whether it is  
better to choose them at random or just use a fixed selection.
The numerical test Sqrt[N[n]]==Round[Sqrt[N[n]] is even faster, but  
it will sometimes produce \false squares\.
Andrzej Kozlowski
===
Subject: Re:  Re: Finding the Number of Pythagorean Triples below a bound
 Hello all,
 emailed
> me.
 To recall, the problem was counting the number of solutions to a
> bivariate polynomial equal to a square,
 Poly(a,b) = c^2
 One form that interested me was the Pythagorean-like equation:
 a^2 + b^2 = c^2 + k
 for {a,b} a positive integer, 0<a<=b, and k any small integer. I was
> wondering about the density of solutions to this since I knew in the
> special case of k=0, let S(N) be the number of primitive solutions
> with
> c < N, then S(N)/N = 1/(2pi) as N -> inf.
 For k a squarefree integer, it is convenient that any solution is
> also
> primitive. I used a simple code that allowed me to find S(10^m) with
> m=1,2,3 for small values of k (for m=4 took my code more than 30  
> mins
> so I aborted it). The data is given below:
 Note: Values are total S(N) for *both* k & -k:
 k = 2
> S(N) = 4, 30, 283
 k = 3
> S(N) = 3, 41, 410
 k = 5
> S(N) = 3, 43, 426
 k = 6
> S(N) = 3, 36, 351
 Question: Does S(N)/N for these also converge? For example, for the
> particular case of k = -6, we have
 S(N) = 2, 20, 202
 which looks suspiciously like the ratio might be converging.
 Anybody know of a code for this that can find m=4,5,6 in a  
> reasonable
> amount of time?

> Yours,
 Titus
> Here is a piece code which utilises the ideas I have described in
>> my previous posts:
>> ls = Prime /@ Range[3, 10];
>> test[n_] :=
>>       Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n],
>> Integers]
>> f[P_, k_] := Sum[If[(w =
>>    a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a,
>>           Floor[Sqrt[P^2 - a^2]]}]
>> g[m_, k_] := f[10^m, k] + f[10^m, -k]
>> We can easily confirm the results of your computations,.e.g. for k=2.
>> Table[g[i,2],{i,3}]
>> {4,30,283}
>> Since you have not revealed the \simple code\ you have used it is
>> hard to tell if the above one is any better. It is however,
>> certainly capable of solving the problem for m=4:
>> g[4,2]//Timing
>> {4779.39 Second,2763}
>> The time it took on my 1 Gigahertz PowerBook was over 70 minutes,
>> which is longer than you thought \reasonable\, so I am still not
>> sure if this is any improvement on what you already have. The time
>> complexity of this algorithm seems somewhat larger than exponential
>> so I would expect that it will take about 6 hours to deal with n=5
>> on my computer, and perhaps 2 weeks to deal with n=6.
>> Andrzej Kozlowski
>
 I mistakenly copied and pasted a wrong (earlier) definition of f.
> Here is the correct one:
 f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
>      1 && (w = a^2 +
>           b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
>            P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]]
 The definition of g is as before:
 g[m_, k_] := f[10^m, k] + f[10^m, -k]
 Andrzej Kozlowski
>
Below is a faster version of the above code. (It owes a significant  
improvement to Daniel Lichtblau, which I borrowed from him without  
his knowledge ;-))
test[n_] :=
     JacobiSymbol[n,
   Prime[Random[Integer, {2, 20}]]] \[CapitalEth] -1 && Element[Sqrt[n],
Integers]
f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
     1 && (w = a^2 +
          b^2 - k) < P^2 && test[w], 1, 0], {a, 1,
           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2+k]]}]]
g[m_, k_] := f[10^m, k] + f[10^m, -k]
The improvement is in the upper bound on a in the sum. Since a is the  
smaller of the two squares whose sum is equal to P^2+k it can't be  
larger than  Floor[Sqrt[(P^2+k)/2]].
Note that you can improve the performance by loosing some accuracy if  
you use a cruder test for a perfect square:
test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w]
    ]
f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] ==
     1 && (w = a^2 +
          b^2 - k) < P^2 && test1[w], 1, 0], {a, 1,
           Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2+k]]}]]
g1[m_, k_] := f1[10^m, k] + f1[10^m, -k]
Let's compare the two cases.
In[7]:=
g[3,2]//Timing
Out[7]=
{89.554 Second,283}
In[8]:=
g1[3,2]//Timing
Out[8]=
{37.376 Second,283}
So we see that we get the same answer and the second approach is  
considerably faster. However:
In[9]:=
g[1,6]//Timing
Out[9]=
{0.008863 Second,3}
In[10]:=
g1[1,6]//Timing
Out[10]=
{0.005429 Second,5}
The correct answer is 3 (returned by the first method). The faster  
method found two false solutions. This should not matter if you are  
interested only in approximate answers (as you seem to be) but it is  
worth keeping in mind.
Andrzej Kozlowski
===
Subject: Re: Finding the Number of Pythagorean Triples below a bound
Hello Andrzej,
The code I used was this:
c1=Sqrt[N[a^2+b^2-k]];
countTriples[m_]:=Block[{a,b,i=0},For[a=1,a<(10^m/Sqrt[2]),a++,
For[b=a,b<10^m,b++,If[c1==Round[c1]&&c1<10^m,i++]]];i]
though m=3 takes about 50 sec (in an old comp) while m=4 will prolly be
around 5000 sec (83 mins). 
-Titus