Carl Love

## 16479 Reputation

7 years, 180 days
Mt Laurel, New Jersey, United States
My name was formerly Carl Devore. I was in the PhD math program at University of Delaware until 2005. I was very active in the Maple community at that time.

## 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.

## Let's correct the original formula first...

Let's correct the original formula to a proper PDF before getting carried away with enthusiasm.

restart;
with(Statistics):
Density := t-> exp(-t^2)/sqrt(Pi):
Check that it's a PDF:
int(Density(t), t= -infinity..infinity);
1
Horizontal shift won't change the integral. But we shift the PDF, not the CDF.
Z:= Distribution(PDF= unapply(Density(t-1), t)):
Now the other aspects can be computed from this automatically by the Statistics package.
CDF(Z,t);
1   1
- + - erf(t - 1)
2   2
Mean(Z);
1

## Let's correct the original formula first...

Let's correct the original formula to a proper PDF before getting carried away with enthusiasm.

restart;
with(Statistics):
Density := t-> exp(-t^2)/sqrt(Pi):
Check that it's a PDF:
int(Density(t), t= -infinity..infinity);
1
Horizontal shift won't change the integral. But we shift the PDF, not the CDF.
Z:= Distribution(PDF= unapply(Density(t-1), t)):
Now the other aspects can be computed from this automatically by the Statistics package.
CDF(Z,t);
1   1
- + - erf(t - 1)
2   2
Mean(Z);
1

## @williamov Use unapply, not Unapply...

Use unapply, not Unapply. It is meaningless with a capital U. Change the line

RealDist := Distribution(Dist(t-1));

to

RealDist := Distribution(CDF= unapply(Dist(t-1), t));

Of course, Markiyan's objection about the original not being a PDF applies here also. I am just trying to offer basic corrections to syntax, not claiming that the above is mathematically valid.

## @williamov Use unapply, not Unapply...

Use unapply, not Unapply. It is meaningless with a capital U. Change the line

RealDist := Distribution(Dist(t-1));

to

RealDist := Distribution(CDF= unapply(Dist(t-1), t));

Of course, Markiyan's objection about the original not being a PDF applies here also. I am just trying to offer basic corrections to syntax, not claiming that the above is mathematically valid.

## Editor dropped characters?...

Preben,

I am having some trouble reading your answer. I wonder if the editor dropped some angle brackets.

I see

SecTer:= CrossProduct(,)

____________________^

and

ProjX := diff(x(t), t, t) = -q*DotProduct(SecTer,,conjugate=false)/m;

_______________________________________^^

and the same problem in the line after that (I tried to used uparrows to indicate the problematic commas.)

## Editor dropped characters?...

Preben,

I am having some trouble reading your answer. I wonder if the editor dropped some angle brackets.

I see

SecTer:= CrossProduct(,)

____________________^

and

ProjX := diff(x(t), t, t) = -q*DotProduct(SecTer,,conjugate=false)/m;

_______________________________________^^

and the same problem in the line after that (I tried to used uparrows to indicate the problematic commas.)

## Precise definition of independent equati...

@Markiyan Hirnyk

It is easier to define independence precisely. Assume that the equations have all been put in the form where the right sides are 0. They are independent if there exists an assignment of numerical values to all of the unknowns such that

1. none of the evaluated left sides are 0
2. none of the evaluated left sides are equal to any other evaluated left side.

PS: Upon sleeping on it, I woke up realizing the above definition is not adequate. I withdraw it, I'll need to think about it some more.

## Precise definition of independent equati...

@Markiyan Hirnyk

It is easier to define independence precisely. Assume that the equations have all been put in the form where the right sides are 0. They are independent if there exists an assignment of numerical values to all of the unknowns such that

1. none of the evaluated left sides are 0
2. none of the evaluated left sides are equal to any other evaluated left side.

PS: Upon sleeping on it, I woke up realizing the above definition is not adequate. I withdraw it, I'll need to think about it some more.

## Maple chose what to eliminate....

Maple chose two of the five unknowns to eliminate. In this case, it is almost obvious that it will pick d and s, but we can't rely on that in general.

﻿