Showing posts with label Mathematica. Show all posts
Showing posts with label Mathematica. Show all posts

Thursday, 13 November 2025

Claude to the Rescue

The number 27983, associated with my diurnal age today, has a property than confers it membership in OEIS A335752: the lexicographically earliest sequence of primes, each two of which sum to a triprime. The sequence begins: 3, 5, 47, 239, 3023, 27983, 53231, ... so let's confirm that these numbers satisfy using 3 as a starting point:
  • \(3 + 5 = 8 =2^3\)
  • \(3 + 47 = 50 = 2 \times 5^2\)
  • \(3 + 239 = 242 = 2 \times 11^2 \)
  • \(3 + 3023 = 3026 = 2 \times 17 \times 89 \)
  • \(3+27983 = 27986 = 2 \times 7 \times 1999\)


The sums above are all triprimes and we will find that any number in the list paired with any other number will produce a triprime. The code for this was included in the OEIS comments but was in Mathematica. I got Gemini to convert it to Python code and then Grok but in both cases the primeomega function was misunderstood with both claiming that primeomega(12) =2 and not 3. 

I didn't fully understand the Python code and was looking to Gemini and then Grok to explain it to me. It was only Claude that succeeded to my satisfaction. Here is the code followed by Claude's explanation. I think Claude will be my go to resource for anything to do with coding.

from sympy import nextprime, primeomega
ps = [3]
for _ in range(10):
    p = nextprime(ps[-1]) 
    while not all(primeomega(x + p) == 3 for x in ps):
        p = nextprime(p)
    ps.append(p)
print(ps)

I'd be happy to explain this code line by line! This code generates a sequence of prime numbers with a special property.

Code Breakdown

from sympy import nextprime, primeomega

Imports two functions from SymPy:
  • nextprime(\(n\)): Returns the next prime number greater than \(n\)
  • primeomega(\(n\)): Returns the total number of prime factors of \(n\) (counting multiplicity). For example, primeomega(12) = 3 because 12 = 2² × 3 (two 2's and one 3)
ps = [3]

Initializes the sequence with the first prime number 3. This list will store our special sequence of primes.

for _ in range(10):

Loops 10 times to find 10 more primes (so we'll have 11 total including the initial 3). The underscore _ is used when we don't need the loop variable.

p = nextprime(ps[-1])

Gets the next prime after the last element in our list. 
  • ps[-1] accesses the last element of the list
  • Initially, this gives nextprime(3) = 5
  • This is our candidate prime to add to the sequence
while not all(primeomega(x + p) == 3 for x in ps):
    p = nextprime(p)

This is the key part! It checks if the candidate prime \(p\) satisfies a special condition:
  • For every number \(x\) already in our list ps
  • The sum \(x + p\) must have exactly 3 prime factors (counting multiplicity)
  • If this condition fails for any \(x\), we move to the next prime candidate
  • This continues until we find a prime where all sums have exactly 3 prime factors
ps.append(p)

Adds the qualified prime to our sequence once it passes the test.

print(ps)

Prints the final sequence** of primes.

Example Walkthrough

Let's trace the first iteration:
1. Start with ps = [3]
2. Try p = 5: Check if 3 + 5 = 8 = 2³ has 3 prime factors ✓ (yes, three 2's)
3. Add 5 to list: ps = [3, 5]
4. Try p = 7: Check 3 + 7 = 10 = 2 × 5 (only 2 factors) ✗
5. Try p = 11: Check 3 + 11 = 14 = 2 × 7 (only 2 factors) ✗
6. Continue until finding a prime that works with both 3 and 5...

This creates a sequence where each new prime maintains the property that when added to any previous prime in the sequence, the result always has exactly 3 prime factors!

*************************

Now if we change the condition to primeomega = 2, we get the following sequence: 3, 7, 19, 139, 859, 8179, 173059, ... where it can be seen there is a big jump from 8179 to 173059. Testing this out using 3 again we see that:
  • \(3+7=10=2 \times 5\)
  • \(3+19=22=2 \times 11\)
  • \(3+139=142=2 \times 71\)
  • \(3+859=862=2 \times 431\)
  • \(3+8179=8182=2 \times 4091\)
  • \(3+173059=173062=2 \times 86531\)
The sequence is included in the OEIS as A113875: slowest growing sequence of primes having the prime-pairwise-average property: if i<j, (a(i)+a(j))/2 is prime. This makes sense in that 2 seems to be the common factor is all of the sums. 

The sequence with additional (very large) terms looks like this:
3, 7, 19, 139, 859, 8179, 173059, 1026199, 1827139, 15828679, 13187242759, 18732483199, 912492556939, 9130567625119

When we set primeomega = 4 we get the following sequence: 3, 13, 137, 1193, 4373, 26017, 33889, 495433, ... which is NOT listed in the OEIS. Again let's use 3 as our test number:
  • \(3+13=16=2^4\)
  • \(3+137=140=2^2 \times 5 \times 7\)
  • \(3+1193=1196=2^2 \times 13 \times 23\)
  • \(3+4373=4376=2^3 \times 547\)
  • \(3+26017=26020=2^2 \times 5 \times 1301\)
  • \(3+33889=33892=2^2 \times 37 \times 229\)
  • \(3+495433=495436=2^2 \times 29 \times 4271\)
That's probably a good place to stop.

Monday, 4 August 2025

Counting Multiplicities in Partitions

I used a new LLM (z.ai) to help me convert Mathematica code to SageMath code and then explain how it all worked. The Mathematica code arose from one of the properties of the number associated with my diurnal age today: 27882. This number qualifies for membership in OEIS A184641.


 A184641: number of partitions of \(n\) having no parts with multiplicity 6.

The sequence begins as follows (permalink):

1, 1, 2, 3, 5, 7, 10, 15, 21, 29, 40, 54, 72, 96, 127, 166, 216, 279, 358, 457, 580, 735, 924, 1159, 1446, 1799, 2228, 2752, 3388, 4158, 5087, 6207, 7551, 9165, 11093, 13401, 16144, 19412, 23286, 27882, 33310, 39727, 47289, 56191, 66647, 78923, 93299

The Mathematica code is:

b[n_, i_] := b[n, i] = If[n == 0, {1, 0}, If[i < 1, {0, 0},
Sum[Function[l, If[j == 6, {l[[1]], l[[1]]}, l]][b[n - i*j, i - 1]], {j, 0, n/i}]]];
a[n_] := b[n, n][[1]] - b[n, n][[2]];
Table[a[n], {n, 0, 50}]

This is incomprehensible to me and so I got z.ai to change this code to SageMath code. I could still not understand the code generated and so I asked the following question:
I really don't understand this code. It doesn't make use of SageMath's Partition function which would seem necessary. Why is the use of this function avoided and is it possible to create more understandable code using this function. I really need to have this code explained step by step. My approach, which is admittedly lengthy, would involve using the Partition function to generate all the partitions for a particular value of \(n\) between 0 and 50. For each value of \(n\), I'd count how many partitions do not contain a part with multiplicity 6 and then move on to the next value of \(n\). Is this approach wrong-headed and too processor intensive. If so, then how can your approach be summarised in layman's terms (like I've done when describing my approach).
Here was the response in part:

This approach is not wrong-headed, but it could be processor-intensive for larger \(n\) because the number of partitions grows exponentially with \(n\). Here's a more understandable implementation using your approach:
def a(n):
    """
    Count the number of partitions of n where no part appears exactly 6 times.
    """
    if n == 0:
        return 1  # The empty partition
    
    count = 0
    # Generate all partitions of n
    for p in Partitions(n):
        # Get the multiplicities of each part
        multiplicities = p.to_exp()
        # Check if 6 is among the multiplicities
        if 6 not in multiplicities:
            count += 1
    return count

# Generate list of a(n) for n from 0 to 50
result = [a(n) for n in range(0, 51)]
print(result)

I then asked:

Can you explain the line "multiplicities = p.to_exp()".

The response was: 

The line multiplicities = p.to_exp() is a key part of the SageMath implementation that deserves a detailed explanation. Let me break it down step by step:

What p.to_exp() Does

    1. Partition Representation:

      • In SageMath, a partition p of an integer n is typically represented as a list of integers in non-increasing order (e.g., [5, 3, 2] for a partition of 10)
      • The .to_exp() method converts this into "exponential form"
    2. Exponential Form:

      • Instead of listing the parts, it lists the multiplicities (how many times each part appears)
      • The output is a list where:
        • The index position represents the part size (starting from 1)
        • The value at that position represents how many times that part appears

Example Walkthrough

Let's take a concrete example to understand this:

Partition: [4, 4, 3, 1, 1, 1] (which sums to 14)

Step-by-step conversion:

    1. Identify all parts and their counts:

      • Part 1 appears 3 times
      • Part 2 appears 0 times
      • Part 3 appears 1 time
      • Part 4 appears 2 times
      • Parts 5 and above appear 0 times
    2. Create the exponential form list:

      • Index 0 (part 1): 3
      • Index 1 (part 2): 0
      • Index 2 (part 3): 1
      • Index 3 (part 4): 2
      • Index 4+ (parts 5+): 0 (but these trailing zeros are omitted)
    3. Result: [3, 0, 1, 2]

Visual Representation

Partition: [4, 4, 3, 1, 1, 1]
Part sizes: 1, 2, 3, 4
Multiplicities:
Part 1: 3 times
Part 2: 0 times
Part 3: 1 time
Part 4: 2 times

Exponential form: [3, 0, 1, 2]

Why This is Useful for Our Problem

In our code, we want to check if any part appears exactly 6 times. With exponential form:

  1. We get a list of all multiplicities
  2. We can simply check if the number 6 appears in this list
  3. If 6 is in the list, it means some part appears exactly 6 times

More Examples

Example 1: Partition with a part appearing 6 times

  • Partition: [6, 6, 6, 6, 6, 6] (six 6's)
  • Exponential form: [0, 0, 0, 0, 0, 6]
  • 6 in [0, 0, 0, 0, 0, 6] → True (we exclude this partition)

Example 2: Partition without any part appearing 6 times

  • Partition: [5, 5, 5, 3, 3, 1]
  • Exponential form: [1, 0, 2, 0, 3]
  • 6 in [1, 0, 2, 0, 3] → False (we include this partition)

Example 3: Partition with multiple parts having different multiplicities

  • Partition: [7, 7, 4, 4, 4, 4, 2, 2, 2, 2, 2]
  • Exponential form: [0, 5, 0, 4, 0, 0, 2]
  • 6 in [0, 5, 0, 4, 0, 0, 2] → False (we include this partition)

Key Advantages of Using to_exp()

  1. Efficiency:

    • Avoids manually counting multiplicities for each part
    • Built-in SageMath function optimized for this operation
  2. Clarity:

    • Directly gives us the information we need (multiplicities)
    • Makes the check for multiplicity 6 straightforward
  3. Completeness:

    • Includes all part sizes from 1 up to the largest part
    • Even parts that don't appear (multiplicity 0) are represented

I then tried the to_exp() out and couldn't get it to work so I asked:

Why is it that when I run the SageMath code in SageMathCell: 
P=[7, 7, 5, 5, 5, 4, 3, 3, 3, 2, 1, 1, 1] print(P.to_exp()) 

I get the error message "AttributeError: 'list' object has no attribute 'to_exp'" 

The response was quite enlightening:

The error occurs because you're trying to call the to_exp() method on a Python list, but to_exp()is a method specific to SageMath's Partition objects. Here's the explanation and solution:

Why the Error Occurs

  1. Python List vs. SageMath Partition Object:
    • When you write P = [7,7,5,5,5,4,3,3,3,2,1,1,1], you're creating a standard Python list
    • Python lists don't have a to_exp() method - this method only exists for SageMath's Partition objects
    • The error message 'list' object has no attribute 'to_exp' tells you exactly this

How to Fix It

You need to convert your list into a SageMath Partition object first:

sage
# Create a Python list
P_list = [7,7,5,5,5,4,3,3,3,2,1,1,1]

# Convert to a SageMath Partition object
P_partition = Partition(P_list)

# Now you can call to_exp()
print(P_partition.to_exp())

Complete Working Example

sage
# Define the partition as a list
partition_list = [7,7,5,5,5,4,3,3,3,2,1,1,1]

# Convert to SageMath Partition object
p = Partition(partition_list)

# Now call to_exp() on the Partition object
multiplicities = p.to_exp()
print(multiplicities)

Expected Output

For the partition [7,7,5,5,5,4,3,3,3,2,1,1,1], the output will be:

[3, 1, 3, 1, 3, 0, 2]

Explanation of the Output

The output [3, 1, 3, 1, 3, 0, 2] represents:

  • Index 0 (part 1): appears 3 times
  • Index 1 (part 2): appears 1 time
  • Index 2 (part 3): appears 3 times
  • Index 3 (part 4): appears 1 time
  • Index 4 (part 5): appears 3 times
  • Index 5 (part 6): appears 0 times (not in the partition)
  • Index 6 (part 7): appears 2 times

Key Points to Remember

  1. Partition Objects vs. Lists:

    • SageMath has special Partition objects with extra functionality
    • Python lists are just basic data structures without these specialized methods
  2. Conversion is Simple:

    • Always convert your list to a Partition object first: Partition(your_list)
    • Then you can use all the SageMath partition methods like to_exp()
  3. Partition Requirements:

    • The list should represent a valid partition (non-increasing integers)
    • SageMath will automatically sort it if needed, but it's good practice to provide it in sorted order

Corrected Version of the Original Code

Here's how to fix the original code snippet you were trying to run:

sage
# Create the partition as a list
P = [7,7,5,5,5,4,3,3,3,2,1,1,1]

# Convert to Partition object
partition = Partition(P)

# Now to_exp() will work
print(partition.to_exp())

This will correctly output [3, 1, 3, 1, 3, 0, 2] without any errors. 

This whole exercise was very useful because it got me using this multiplicity function that I was completely unaware of. It's important to note that the Python list being converted to a SageMath partition is put into descending order. For example, consider the list P = [5, 6, 5, 4, 4, 7, 4, 3, 1, 1, 2, 3, 1]. Using "sorted" before the list we NOT work because the elements of the list are put in ascending order. What's needed is "reversed(sorted)" as shown below:

P=[5, 6, 5, 4, 4, 7, 4, 3, 1, 1, 2, 3, 1]
Q=Partition(reversed(sorted(P)))
R=Q.to_exp()
print(Q,"-->",R)

[7, 6, 5, 5, 4, 4, 4, 3, 3, 2, 1, 1, 1] --> [3, 1, 2, 3, 2, 1, 1]

However, even though the partition members must be arranged in descending, the output of the multiplicities is in ascending order. Thus we have three 1's, one 2, two 3's and so on. This can be a little confusing.

Here is another application of the to_exp() function, this time used to find the numbers of partitions of \(n\) that contain no multiplicities (OEIS A000009):

# Partitions of n that have no multiplicity higher than 1

L=[]
for n in [0..40]:
    P=Partitions(n)
    count=0
    for p in P:
        multi=p.to_exp()
        # Check if all multiplicities are <= 1
        if all(m <= 1 for m in multi):
            count+=1
    L.append(count)
print(L)

[1, 1, 1, 2, 2, 3, 4, 5, 6, 8, 10, 12, 15, 18, 22, 27, 32, 38, 46, 54, 64, 76, 89, 104, 122, 142, 165, 192, 222, 256, 296, 340, 390, 448, 512, 585, 668, 760, 864, 982, 1113] 

I'm thankful to z.ai suggestion of "if all(m <= 1 in multi)" because it's more succinct than my normal approach which is to set a variable OK to 1 and then reset it to 0 if a condition is not met. After that I test if OK is still equal to 1 and, if so, then I increment the count. This single line replaces the seven lines previously required:

OK=1
for m in multi:
    if m >1:
        OK=0 
        break 
if OK==1:
    count+=1

I need to remember that function for future use because my coding is ... well, terrible.