Carl Love

## 25249 Reputation

10 years, 230 days
Himself
Natick, Massachusetts, United States
My name was formerly Carl Devore.

## Please clarify "power is limited"...

When you say "their power is limited ... [to], for instance, no more than 3", do you mean that the total degree of the polynomials is at most 3? Or do you mean that the degree, or power, of each individual variable is at most 3?

## What do you mean by "solve"?...

What do you mean by "solve"? Do you mean "approximate"? What is the "it" about whose truth you ask?

You've given some Maple commands, which work correctly. But what's the point of only using three terms of your series? Why aren't you satisfied with the output?

## No, it's not too large...

To address directly the OP's issues:

mike_a wrote:

> I noticed some strange things with the list/table carrying my complex points. If I ask for individual values, like T, it won't give the correct value, but if I ask for a range, like T[1..20], it will display the correct values.

You'll have to post some code so we can figure out what's going on. But I'll tell you right off that your problem has nothing to do with the size of your data set.

> I get an error when trying to use complexplot saying either invalid range, if I don't have a range for the 2nd and 3rd argument, or else Error, (in plot) procedure expected, as range contains no plotting variable. I assume that there's some problem with the table being too large and not a problem with the complexplot command

No, it has nothing to do with the table being too large. But you'll need to post some code so that we can figure out the actual problem.

First work out these indexing problems on a much smaller version of your plot before moving on to the million-point plot.

## Summary...

To summarize what you said, for those who may have trouble reading it: You are trying to solve, using a least-squares method, a linear system A.X = B, where A has dimensions 1201 x 800 and B is a column vector. You have a tried an ad hoc method that involves explicitly inverting an 800 x 800 matrix, but ran out of memory.

You wrote:

> A^T*A*X=A^T*X.  let M=A^T*A.  M*X=A^T*X  finally  X=M^-1 * A^T*X

That's supposed to be B, not X, on the right sides of those equations! I hope that was just a typo in your post and not in your Maple code!

Several questions before we go any further:

1. Is A a hardware floating point matrix?
2. What percentage of the elements of A are non-zero? A rough estimate is fine.
3. Have you tried LinearAlgebra:-LeastSquares?
4. Can you make a smaller version of this problem that we can practice with before doing the full 1201 x 800?

## Plot of Preben's one-term approximation...

A plot shows that Preben's one-term approximation is excellent even for x close to -2, and the original is visually indistinguishable from the approximation for x < -3.2. The plot below uses A and res as defined in Preben's code. If Digits, numpoints, and adaptive are not reduced from their default values (as I have done below), Maple will take a very long time to do the numerical integrations.

Digits:= 6:
plot([op(1,res), A], x= -3.3..-2.1, numpoints= 50, adaptive= 4); I had originally plotted down to x = -10, but then it difficult to even see that there are two curves because they are so close.

## Plot of Preben's one-term approximation...

A plot shows that Preben's one-term approximation is excellent even for x close to -2, and the original is visually indistinguishable from the approximation for x < -3.2. The plot below uses A and res as defined in Preben's code. If Digits, numpoints, and adaptive are not reduced from their default values (as I have done below), Maple will take a very long time to do the numerical integrations.

Digits:= 6:
plot([op(1,res), A], x= -3.3..-2.1, numpoints= 50, adaptive= 4); I had originally plotted down to x = -10, but then it difficult to even see that there are two curves because they are so close.

## Are you sure of that formula?...

Are you sure that you didn't switch n and m in your big-O formula? Also, do you mean to compute separately the GCDs of m pairs of numbers? or to compute the overall GCD of a set of m numbers?

## Iterative versus sequential output and b...

I think the key to the timing differences w.r.t. production of output might be in the line

K[cnt]:= seq([j[1..-2][], Div[k], Div[-k]], k= `if`(i=2, 1, BinarySearch(j[-2]))..Nsq)

where all the splits of [..., n] are generated in a single call to seq. If the output was being generated one at a time, as with an iterator, I couldn't use that efficient seq. In the Iterator package, could you generate iterates in bursts like this and store them in a buffer? I know one goal of the package is minimal memory usage, but in this case the seq is generating only a small portion of the overall output.

## With binary search and local remember ta...

Before posting my last version, I had tried with ListTools:-BinarySearch, but that consumed about half of the time. This time, I implemented a binary search tailored to this algorithm. As Joe said, it's probably not worth it; but it's probably not a significant waste either.

A much more significant change is that I made it call numtheory:-divisors only once and sort only once. All subsequent sorted sublists of divisors are made my filtering the this one list. And the filtered lists are stored in a local remember table.

It seems that there is considerable interest in this algorithm. It certainly fascinates me. It amazes me that simply keeping the lists of factorings carefully sorted and only splitting the last element is sufficient to generate all the factorings. Joe: Is your method for generating all partitions of a multiset similar to this?

Factorings:= proc(n::posint, m::posint:= 0)
local
L, K, Div, i, j, k, Nsq, cnt
,Omega:= numtheory:-bigomega(n)  # max # of factors
,DIVS:= sort([(numtheory:-divisors(n) minus {1,n})[]])  # not changed during this proc's run.

# Get sorted list of divisors of divisors by filtering DIVS; also returns pointer to midpt of list.
,divisors:= proc(n::posint)
option remember;
local D:= select(x-> irem(n,x)=0, DIVS)[1..-2];  # remove n itself from end.
D, iquo(nops(D)+1, 2)  # Extra 1 is for square root.
end proc

# Find first index k into Div[1..Nsq] such that Div[k] >= x.
,BinarySearch:= proc(x)
local lo, hi, m;
if Nsq=0 then  return 1  elif Div[Nsq] < x then  return Nsq+1  fi;  # fall-thru cases.
lo:= 1; hi:= Nsq;
do
m:= iquo(hi+lo, 2);
if Div[m] = x then  return m  elif x < Div[m] then  hi:= m-1  else  lo:= m+1  fi;
if lo > hi then  return lo  fi
end do
end proc
;
if m > Omega then  return []  elif m > 0 then  Omega:= m  fi;

L:= [[n]];  # Initialize loop w/ sole factoring of length 1.
# Avoid wasted effort of filtering divisor list on the first (i=2) pass.
divisors(n):= DIVS, iquo(nops(DIVS)+1, 2);

for i from 2 to Omega do  # i is number of factors.
cnt:= 0;  # cnt is index for list-building table.
# j is a single factoring for previous i, so j is a list.
for j in L[i-1] do
cnt:= cnt+1;
# Split last, and greatest, member of list, if not prime (prime case falls thru).
Div, Nsq:= divisors(j[-1]);
# j[-1] = Div[k]*Div[-k] for any k, since Div is sorted.
# Only use splits that maintain the list in order.
K[cnt]:= seq([j[1..-2][], Div[k], Div[-k]], k= `if`(i=2, 1, BinarySearch(j[-2]))..Nsq)
end do;
# Convert table K to list L[i].
L[i]:= [seq](K[j], j= 1..cnt)
end do;

# Convert table L to list for final output.
`if`(m = 0, [seq](L[i][], i= 1..Omega), L[m])
end proc;

With adaptive set to false, the plot lacks the rich detail I want unless numpoints is set to 2^15 or greater. This is a single closed curve plotted in a single color, so all the action is on the boundary (the whole thing is its boundary). I'm willing to cut back to numpoints= 2^12, adaptive= 4. A little time can be saved, without significantly impacting quality, by truncating the series at 9 terms (N=8) rather than the 17 terms I used above.

Go ahead and use it as an example, and give me credit for the X and Y functions. You could use 2^15 evenly spaced points and do it very quickly in evalhf.

Here is a much larger jpeg showing more detail. This is with N=8, numpoints= 2^12, adaptive= 4. Blow it up and see how the inner curlicues look like flowers. ## More modifications...

It's a great algorithm.  I've made some modifications to the code in addition to Joe's, keeping the basic algorithm the same. You may notice another factor-of-two time improvement on some large factorings such as all 92,213 factorings of our friend 711*100^3.

1. I eliminated the repetitive deep indexing by letting j become a factoring itself rather than an index into the list L[i-1] of factorings.
2. I rolled the factorings of length two case into the main loop.
3. I eliminated all divisions by always sorting the list of divisors. That is, if D:= sort([divisors(N)[]]), then for all k, we have n = D[k]*D[-k].
4. Because of (3), I only scan the first half of each list of divisors, plus one more in case of a perfect square (we know it's a square if the number of divisors is odd, so I don't need to check).

Factorings:= proc(n::posint, m::posint:= 0)
uses NT= numtheory;
local
L, T, Div, i, j, k, Nsq, idx
,Omega:= NT:-bigomega(n)
;
if m > Omega then  return []  elif m > 0 then  Omega:= m  fi;

L:= [[n]];

for i from 2 to Omega do  # i is number of factors.
idx:= 0;
# j is one factoring for previous i, so j is a list.
for j in L[i-1] do
# Split last, and greatest, member of list, if not prime.
Div:= sort([(NT:-divisors(j[-1]) minus {1,j[-1]})[]]);
Nsq:= iquo(nops(Div)+1, 2);  # Extra 1 is for square root.
# Only use splits that maintain the list in order.
if i=2 then  k:= 1  else  for k to Nsq while Div[k] < j[-2] do od  fi;
idx:= idx+1;
# j[-1] = Div[k]*Div[-k] for any k, since Div is sorted.
T[idx]:= seq([j[1..-2][], Div[k], Div[-k]], k= k..Nsq)
end do;
# Convert table T to list L[i].
L[i]:= [seq](T[j], j= 1..idx)
end do;

# Convert table L to list for final output.
`if`(m = 0, [seq](L[i][], i= 1..Omega), L[m])
end proc:

## symbols won't expand...

Kitonum: The reason to use ``(...) rather than convert(..., symbol) is that expand will remove the ``.

## symbols won't expand...

Kitonum: The reason to use ``(...) rather than convert(..., symbol) is that expand will remove the ``.

## Post example worksheet...

Yes, it would be nice and would make sense if simplify were idempotent.