Carl Love

Carl Love

26513 Reputation

25 Badges

11 years, 190 days
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are replies submitted by Carl Love


Here is a procedure to reconstruct the original from the output of PolynomialReduce(P,R,V). It only requires one preliminary step to your expand(inner(...)) + ..., and it needs to be passed the elimination order that was used.

PolyRecon:= proc(Q::list(polynom), Rm::polynom, Z::list(polynom), V::list(name))
local k, z:= <Z>, v;
    for k,v in V do z[k+1..]:= rem~(z[k+1..], z[k], v) od;
    expand(<Q>^+.z + Rm)
end proc
(P,R,V):= (x^5 - x*y^6 - x*y, [x^2 + y, x^2 - y^3], [y, x]):
[    5    3  4    5  3    7  2    9     / 10    \      7    3]   3
[-x y  + x  y  - x  y  + x  y  - x  y + \x   - 1/ x, -x  + x ], x 

PolyRecon(%, R, V);
                            6    5      
                        -x y  + x  - x y
V:= [x,y]:
                  [ 3           3      ]      
                  [x  - x y, x y  - x y], -x y

PolyRecon(%, R, V);
                            6    5      
                        -x y  + x  - x y


@sursumCorda I wouldn't expect your reconstruction process to return the original, just something equivalent to the original under the given reductions. It seems to me that it does that. In other words, for any P, R, V such that PolynomialReduce(P,R,V) returns a result, I think that it'll always be true that

PolynomialReduce(expand(inner(R, PolynomialReduce(P,R,V)[1])), R, V)[2] = 0

Note that for each iteration of my loop, all future divisors are also reduced.

@sursumCorda The elimination order matters a great deal in these problems. For the case shown, 
rem(y^3 - y, x^3 - x, y) = 0, so for the second loop iteration, the divisor is 0. I suppose one could trap for this error and try a different order in that case.

(My knowledge of this area of math (I guess it's Commutative Ring Theory?) is limited.)

Here is a worksheet that summarizes all of the above.


FamN:= [3,2,2,2,1]:  #reduced family sizes for test problem

nf:= nops(FamN):

PS:= [0, seq[scan= `+`](FamN)];  #their partial sums

[0, 3, 5, 7, 9, 10]

FamS:= (`{}`@`$`@`..`)~(PS[..-2]+~1, PS[2..]);  #partition

[{1, 2, 3}, {4, 5}, {6, 7}, {8, 9}, {10}]

It:= Iterator:-Permute(PS[-1]):

# For efficiency, split the Iterator's output array into aliases for
# the image of each family under a given permutation:
sp:= output(It):
FamA:= [seq](ArrayTools:-Alias(sp, PS[i], [FamN[i]]), i=

# Does the given permutation swap families i and j?
Swapped?:= (i,j)-> local x;
    andseq(x in FamS[j], x= FamA[i]) and andseq(x in FamS[i], x= FamA[j])

SwapJ:= (2,3):  #indices of families to swap

V:= 0:  #counter of permutations with the desired property

CodeTools:-Usage(in It do if Swapped?(SwapJ) then V++ fi od):

memory used=1.04GiB, alloc change=26.27MiB, cpu time=11.08s, real time=10.04s, gc time=1.83s



V/10!;  #probability


# Now, do the same thing with fixed-size set partitions, which are often
# easier to work with because there are far fewer of them than permutations.

CountPartitions:= (S::list(posint))->
    add(S)!/mul([S[], rhs~(Statistics:-Tally(S))[]]!~)

# My flawed first-posted attempt was analogous to this:
CountPartitions(FamN[[({$nf} minus {SwapJ})[]]])
/ CountPartitions(FamN);


# The correction is:
CountPartitions(FamN[[({$nf} minus {SwapJ})[]]])
/ CountPartitions(FamN) / numboccur(FamN, FamN[SwapJ[1]])!


# Now apply what we've learned above to the original problem:
FamN:= [3,4,4,4,1]:
CountPartitions(FamN[[({$nops(FamN)} minus {SwapJ})[]]])
/ CountPartitions(FamN) / numboccur(FamN, FamN[SwapJ[1]])!


# So, the probability that this happens for any one of the 3 possible pairs
# of equal-size families, and happens for that same pair in the next
# drawing is:
p1:= evalf(3*%^2);


# To put that into perspective, the probability of a royal flush in a 5-card
# poker hand with no card exchanges is
p2:= evalf(4/binomial(52,5));






@sursumCorda Here is a procedure that might do what you want. Let me know. I've made it so that the elimination order is an optional third argument:

PolynomialReduce:= proc(
    p::polynom, Z::list(polynom), 
    ord::list(name):= [indets(p, And(name, Not(constant)))[]]
local r:= p, v:= ord, z:= Z, k, d;
    if nops(v) < nops(z) then 
        v:= [v[], (indets(z, And(name, Not(constant))) minus {v[]})[]]
    [for k to nops(z) do
         z:= rem~(z[2..], (d:= z[1]), v[k]);
         quo(r, d, v[k], 'r')
end proc
PolynomialReduce(x^2+y^2, [x-y, y+a]);
                   [x + y, -2 a + 2 y], 2 a 

I withdrawal the above Reply. Indeed, my probabilities were too large by a factor of 6. Although my procedure CountPartitions is correct, one does need to consider the order of the images when computing the probability. I confirmed this by doing the analogous problem with the sizes changed to [3,2,2,2,1] and enumerating all 10! ~ 3.6 million permutations, which only takes 11 seconds.

@mmcdara @Christopher2222 The anomalous extra multiplier of 1/6 (exactly) that you're both getting is due to subtly and unawarely implying that the set images of the families under the permutations are in some family-to-family order. In truth, the permuted images of families BCD can appear in any of 3! = 6 orders.

That is why my counting procedure CountPartitions divides out the factorials of the rhss of the Statistics:-Tally of the block sizes S


@acer was only using posint to refer to the type of the suffix on the variable. For example, in _Z2, the is the suffix, and it happens to be a posint. This has nothing to do with what the variable is "assumed to be", which is almost always integer (negative, zero, or positive) for variables beginning _Z. Other common automatically generated variables returned by solve are prefixed _NN (assumed to be a nonnegative integer) and _B (assumed to be 0 or 1). 

Things that are formally assumed about variables are "properties", not "types". It's very easy to confuse properties and types. They are quite similar, but not exactly the same thing. Further adding to the confusion is that many properties and types have the same names, for example integer.

The basic command that checks whether something has a given property (analogous to the type command for types) is is. For example,

solve({sin(x)=0}, {x}, allsolutions);
{x = Pi _Z2~}
is(_Z2, integer);
Originally _Z2, renamed _Z2~:
  is assumed to be: integer

The about command returns the information also, but in a format that's difficult to use programmatically because it's intended for direct human reading.

Finally, note that objects can satisfy many types simultaneously. The command whattype only returns the most fundamental, and I usually discourage using it. Likewise (this is probably more obvious), variables can simultaneously have multiple properties or assumptions.

Hvae a look at the package SignalProcessing (help page ?SignalProcessing). It might have some commands to help with that. I don't know enough about that branch of math to say for sure.


@acer was only using posint to refer to the type of the suffix on the variable. For example, in _Z2, the is the suffix, and it happens to be a posint.

I completely rewrote the rest of this Reply to emphasize and elaborate upon the type/property distinction, and it's further down in this subthread now.

@mmcdara You wrote:

  • Maybe I made a wrong interpretation of the term "exchange" while considering that if someone, X, recieves the gift that someone else, Y, bought, then Y recieves X's.

I think that you are interpreting that part incorrectly. The desired permutations are not necessarily symmetric (i.e., of order 2), as can be seen from the final output of @Christopher2222 's most-recent Reply.

@Christopher2222 I computed the exact probability of families B and C swapping in one year, and in successive years. The probability is of course very small, but it's not quite as small as I was guessing. Unlike the earlier probability that I computed exactly, there's a relatively simple formula for this one. It can even be done with pen and paper, which I did yesterday while I was at a party without access to Maple.

I start with a formula for the number of fixed-size set partitions using a specified list of family sizes (in this case [3,4,4,4,1]).

(* Formula to count the number of fixed-size set partitions with the
specified block sizes. This returns an inert expression (for educational
purposes) that can be evaluated with the value command.
CountPartitions:= (S::list(posint))->
    %factorial(add(S)) %/ 
    mul(%factorial~(subs(1= (), [S[], rhs~(Statistics:-Tally(S))[]])))

So, for the case at hand,
             16! / (3!^2 * 4!^3)

Now let's consider all set partitions where families B and C swap (in exactly the way that you described). The setwise image of B under the corresponding permutations is C, and the setwise image of C is B. So from the point of view of fixed-sized set  partitions, the complete count is the same as the count of the partitions of the remaining 8 people into blocks of sizes [3,4,1].

CountPartitions([3,4,1]);  value(%);
             8! / (3! * 4!)

So  the probability of this happening in any given drawing is
280 / 42042000;
             1 / 150150

and the probability of that happening in two specified drawings is
%^2;  evalf[2](%);
             1 / 22545022500
             4.4 * 10^(-11)

To be fair, we should consider the probability of any pair of families swapping and that same pair swapping in the next drawing. In this case, that's only possible for 3 pairs: (B,C), (C,D), and (B,D). No two of these pairs can swap at the same time, so the probabilities are simply 3 times those computed before, with the final number becoming 1.3 * 10^(-10). I still consider that "infinitesimal"! To put that in perspective, you'd be 11566 times more likely to get the highest possible poker hand---a royal flush---in a 5-card deal with no extra cards than to have the same two families paired 2 years in a row.

@Christopher2222 There's also a command combinat:-randperm, which is easier to apply to this problem than randcomb.


Here I attempt to describe the problem with a degree of mathematical abstraction and less reliance on social and anthropological constructs that might cause confusion of the math.

We have a set of 16 people partitioned into 5 subsets, here called "families". Using some randomization method (colloquially, traditionally, and metaphorically referred to as "drawing names from a hat") we generate a permutation of the 16 so that each person is assigned another from the group to give a gift to. This is not intended to be a symmetric exchange (a transposition), but that is allowed if that's the way the random assignment turns out.

Of course, there is no fun in this if a person is assigned to themself. So we want the permutation to be a derangement[*1] in the standard mathematical sense. To further increase the level of fun, social bonding, etc., we impose a stronger restriction that I'll call a generalized derangement: No person should be assigned to someone from their own partition block (i.e., family).

So, the question is What is the probabilty that a random permutation satisfies this restriction?

I've heard of this tradition being called "Secret Santa".

[*1] This usage of derangement is strictly a mathematical type of permutation (specifically, it's a permutation with no fixed points); it has no connection to the state of mental illness of the same name.

@Christopher2222 I haven't computed it, but off the top of my head I'd guess that the probabilty that 2 of the same-size families would reciprocally exchange given a random drawing is on the order of 10^(-11), so the probability of that happening 2 years in a row is on the order of 10^(-22). As I mentioned in my exact computation, the number of permutations is very close to 2.1×10^(13), aka 21 trillion.

Ah! This is forensic mathematics! 

First 9 10 11 12 13 14 15 Last Page 11 of 689