Applications, Examples and Libraries

Share your work here

I was wondering whether a water molecule exhibits the intermediate axis theorem (IAT) effect.

I was also curious if a tool like MapleSim could be used for simulations, as it is designed with technical applications in mind. (It would have been handy to enter parameters in Angstroms and atomic mass, but MapleSim only provides length in mm and mass in grams for small objects.)  

The atoms of the H2O molecule form an obtuse triangle. An obtuse triangle exhibits the effect when rotated about its intermediate axis (horizontal axis in the image below) provided it consists of equal point masses (more details here) which is not the case for oxygen and hydrogen. 

Building and running a model is quick since geometry and masses are available online. Two molecules are modeled in the attached file with parameters from different sources. The left model in the image below represents the atoms as three discrete point masses without individual rotational inertia. The system's total inertia is derived entirely from the spatial distribution of these masses. On the right, the entire molecule is modeled as a single lumped mass located at the center of mass (CoM), with the molecule’s full rotational inertia tensor applied at that point.

IAT_H2O.msim

Both models are rotating about the assumed intermediate axis at 3 THz (3 trillion rotations per second), which is about two orders of magnitude less than near-infrared light (just beyond the visible spectrum). To generate the IAT effect an orthogonal tiny initial “kick” of 0.001 Hz was added.

Visualizing was not so easy since the models are in the picometer range. It was difficult to locate the molecules by zooming with the mouse wheel (the button “Fit scene” on the 3-D Playback window did not work). Running an animation to visualize the IAT effect was not possible because the display speed could “only” be slowed down by a factor of 2^31 (see A in the image below). However, the slider C to move back and forth in time worked. The time (B) was always rounded to zero due to the very short simulation time span of 1E-10 seconds (i.e. 10 picoseconds).

Export of an animation movie worked better. Over the displayed timespan of 10 ps the molecules flip 2 times back and forth. This corresponds to a flip frequency of 200 GHz.

Time is still not displayed in the movie but that is pretty much all what did not work. Overall, a good performance. Also the numerical results did not show signs of loss of fidelity despite the atomic scale  of the objects and the very short time span.

In reality, the rotation of a dipole results in the emission of electromagnetic radiation, which dissipates energy and decelerates the molecule. How can this damping effect be modeled in MapleSim? Someone has an idea?

For decades, Maple has been built around one of the world’s most powerful mathematics engines—helping students, educators, engineers, and researchers explore ideas, solve complex problems, and communicate mathematics clearly.

Maple 2026 builds on that foundation with major advances in the math engine, expanding the kinds of problems Maple can solve while improving reliability and performance.

At the same time, Maple 2026 introduces new AI-powered tools that help you work faster—finding commands, generating visualizations, explaining concepts, and helping you explore ideas. The key difference is that these tools sit on top of Maple’s math engine, so the results are grounded in real computation rather than guesswork.

If you’ve been following along with our recent Mathy teaser videos and sneak peek posts, you may already have seen hints of some of these features. Now I’m excited to finally share them in full.

One of the most exciting additions in Maple 2026 is the new AI Assistant.

AI tools are incredibly useful for exploring ideas, writing code, and learning new topics. But when the mathematics becomes more involved, relying on AI alone can be risky. The Maple AI Assistant brings those productivity benefits into Maple while keeping the mathematics grounded in Maple’s trusted computation engine.

You can ask the AI Assistant questions in natural language and have it help you:

  • find Maple commands or formulas
  • generate Maple code
  • create visualizations
  • explain mathematical concepts
  • draft examples, worksheets, or reports

Because Maple performs the underlying computations where appropriate, the results are grounded in Maple’s powerful math engine. The AI Assistant becomes a productivity partner that helps you accomplish tasks in Maple faster and more easily, combining the flexibility of AI with mathematics you can trust.

Watch the AI Assistant in action.

 
Turn Documents into Live Mathematics

Another feature I’m particularly excited about is Document Import.

Many of us have years of mathematical content stored in PDFs, lecture notes, journal articles, slides, or even handwritten pages. Traditionally these documents are static—you can read them, but you can’t interact with the mathematics inside them.

With Maple 2026, that changes.

Document Import allows Maple to convert many document formats—including PDFs, DOCX files, and presentations—into Maple worksheets where the mathematics becomes live and executable. 

The image below illustrates the transformation.

On the left (“Before”), scribbled handwritten notes from a Calculus III lecture were saved in a Word document. The notes include hand-drawn sketches, formulas, and written explanations.

After importing the document into Maple (“After”), the mathematical expressions were recognized and converted into live, editable Maple mathematics. The text was preserved, and the hand-drawn sketches were retained as images. The resulting worksheet supports evaluation, editing, and further computation.

Once imported, you can:

  • evaluate expressions
  • modify formulas
  • extend derivations
  • add visualizations
  • explore variations of the mathematics

Instead of recreating examples from scratch, you can bring existing material directly into Maple and start exploring.

While the new AI features are exciting, the heart of Maple has always been its mathematics engine—and Maple 2026 delivers significant advances here.

One particularly notable improvement is Maple’s expanded ability to solve linear recurrence equations. Through improvements to the rsolve command and major extensions to the LREtools package, Maple can now solve dramatically more recurrence relations than before, including many third- and fourth-order cases that were previously beyond reach.

In fact, Maple can now fully solve over 94% of the 55,979 entries in the Online Encyclopedia of Integer Sequences (OEIS) that that can be shown to satisfy a linear recurrence relation. These advances reflect ongoing research into linear difference equations and their algorithmic implementation in Maple, continuing Maple’s long tradition of advancing the state of computer algebra.

Beyond recurrence solving, Maple 2026 includes many improvements across its core symbolic and numeric algorithms. Maple’s assumption system has been strengthened to improve reasoning under mathematical assumptions, and enhancements to the simplify, combine, and evalc commands allow Maple to produce more compact and mathematically natural forms for a wider range of expressions.

There are also improvements to Maple’s differential equation solvers, polynomial system solving, and numerical solving routines such as fsolve, along with updates to other foundational parts of the math library used throughout the system.

Taken together, these improvements expand the range of problems Maple can solve and improve the robustness, correctness, and efficiency of the results.

Maple has always offered extensive control over plotting options, but achieving consistent visual styling across multiple plots could require specifying many settings each time.

Maple 2026 introduces Plotting Themes, which allow you to define a plotting style once and apply it across many plots with a single option.

Themes make it easy to maintain consistent visual styles in worksheets, teaching materials, reports, and publications, while still allowing individual plots to override specific options when needed.

The image below shows an example of creating and applying a custom plotting theme. 

 

Maple continues to be widely used in classrooms around the world, and Maple 2026 includes several improvements designed to support teaching and learning.

The Check My Work system has been enhanced so Maple can recognize a wider variety of valid student solution steps and provide more accurate feedback.

Maple 2026 also improves the generation of similar practice problems, making it easier to create variations of a problem while preserving its mathematical structure.

In addition, Maple’s step-by-step solutions have been expanded to support more types of expressions, helping students better understand the reasoning behind the mathematics they’re learning.

Maple 2026 also introduces improvements for developers building advanced applications, along with performance enhancements across the system.

One particularly interesting addition is the new VectorSearch package, which implements a vector database directly inside Maple.

If you’re not familiar with vector databases, one way to think about them is through recommendation systems like Netflix or Spotify. Each movie or song can be represented by a vector containing thousands of numbers describing its characteristics—things like genre, pacing, or mood. When you watch something, the system finds other items whose vectors are closest to it, which is how recommendations are generated.

With the new VectorSearch package, Maple can store thousands (or more) of vectors and efficiently find the ones most similar to a given vector. This makes it easier to build applications involving machine learning, data analysis, and modern AI workflows directly in Maple.

Maple 2026 also delivers significant performance improvements. For example, operations involving quantities with units have been greatly optimized—some computations now run over 90 times faster, making Maple even more efficient for engineering and scientific workflows.

Maple 2026 also expands the benefits available through the Maplesoft Elite Maintenance Program (EMP). The new benefits include access to additional Maplesoft products and services:

  • Maple Learn, the online environment for teaching and learning mathematics
  • Maple Calculator Premium, bringing the power of Maple to your phone with full access to features like Solution Steps and Check My Work
  • Maple MCP, which allows you to connect Maple’s math engine to external AI tools so they can produce mathematical results you can trust

These additions extend Maple beyond the desktop, giving users powerful tools for learning, teaching, and exploring mathematics across web and mobile platforms, as well as through integrations with external AI tools.

This post only scratches the surface of what’s new in Maple 2026. There are many more improvements across the math library, programming tools, and performance.

To learn more about all the new features and enhancements in Maple 2026, visit the What’s New in Maple page on our website.

 

 

I am very pleased to announce that Volume 6 Number 1 of Maple Transactions has been published.  This is a Special Issue on Matrices and Polynomials in Computer Algebra, and the Guest Editors (our first ever!) were Marc Moreno Maza and Tomas Recio.  There are still some papers that are expected to be added to the issue when they come in, but at this moment there are 8 papers there for you to read (and a description of the issue in the Front Matter section, by the Guest Editors).

A link to this Special Issue

This post was inspired by the following discussion thread  https://mapleprimes.com/questions/242266-Count-The-Number-Of-Paths , which considered the problem of finding all Hamiltonian paths on an integer lattice in R^2 that connect two distinct vertices. The  AllPaths  procedure solves a more general problem: it finds all self-disjoint paths connecting two distinct vertices not only in the plane  R^2  but also in space  R^3 . Of course, it also finds all Hamiltonian paths or allows one to determine their absence. The procedure does not use commands from GraphTheory package (only direct manipulation of sets and lists).
Required parameters of the procedure: is a set or list of lattice vertices specified by their coordinates, Start and Finish are the initial and final vertices. Optional parameter R (defaults it's NULL  if all paths are searched) and any symbol (if only Hamiltonian paths are searched). S can be either a rectangular integer lattice or a union of several such lattices.

Code of the procedure:

restart;
AllPaths:=proc(S::{set(list),listlist},Start::list,Finish::list,R::symbol:=NULL)
local N:=nops(S), S1:=convert(S, set), L, n, m, k, i, j, s, p, q, P, a, b, c;

L:={[Start]};
for n from 2 to N do

if R=NULL then

P:='P';
m:=nops(L);
for k from 1 to m do
if nops(Start)=2 then
i,j:=L[k][-1][];
s:={[i-1,j],[i,j+1],[i+1,j],[i,j-1]} else
a,b,c:=L[k][-1][];
s:={[a-1,b,c],[a,b+1,c],[a+1,b,c],[a,b-1,c],[a,b,c-1],[a,b,c+1]} fi; 
s:=`intersect`(s,S1) minus convert(L[k],set);
if s={} and L[k][-1]=Finish then P[k]:=L[k] else
if s={} and L[k][-1]<>Finish then P[k]:=NULL else
P[k]:=`if`(L[k][-1]=Finish,L[k],seq([L[k][],s[i]],i=1..nops(s)))  fi; fi;
od;
L:=convert(P,set) else

P:='P';
m:=nops(L);
for k from 1 to m do
if nops(Start)=2 then
i,j:=L[k][-1][];
s:={[i-1,j],[i,j+1],[i+1,j],[i,j-1]} else
a,b,c:=L[k][-1][];
s:={[a-1,b,c],[a,b+1,c],[a+1,b,c],[a,b-1,c],[a,b,c-1],[a,b,c+1]} fi;
s:=`intersect`(s,S1) minus convert(L[k],set);
if n<N then 
if s={} or L[k][-1]=Finish then P[k]:=NULL else
P[k]:=seq([L[k][],s[i]],i=1..nops(s)); fi else 
if L[k][-1]=Finish and s<>{} then P[k]:=NULL else P[k]:=seq([L[k][],s[i]],i=1..nops(s));
fi; fi;  
od;
L:=convert(P,set)

fi; od;

L;
end proc:


Examples of use.

In the first example from the post above, we find the number of Hamiltonian paths in 

L:=CodeTools:-Usage(AllPaths({seq(seq([i,j], i=1..11), j=1..3)}, [2,2], [10,2], 'H')):
nops(L);

   

In this same example, we find the number of all paths from A to B and the possible lengths of these paths.

L:=CodeTools:-Usage(AllPaths({seq(seq([i,j], i=1..11), j=1..3)}, [2,2], [10,2])):
nops(L);
map(t->nops(t), L);
L1:=select(t->nops(t)=%[-1], L):
nops(L1);

  

In the following example, we find the number of all paths, as well as the number of Hamiltonian paths, and animate these paths (total 24 one's).

S:={seq(seq([i,j],i=1..5),j=1..3)} union {seq(seq([i,j],i=4..7),j=3..5)}: A:=[1,1]: B:=[7,5]:
P:=plots:-display(plots:-pointplot(S, symbol=solidcircle, color=blue, symbolsize=15, view=[0..7.5,0..6.5], size=[600,500], scaling=constrained), plots:-textplot([[A[],"A"],[B[],"B"]], font=[times,bold,22], align=[left,above])):
L:=AllPaths(S,A,B):
nops(L);
map(t->nops(t), L);
L1:=select(t->nops(t)=%[-1], L):
nops(L1);
plots:-animate((plots:-display)@(plottools:-curve),[L1[round(a)], color=red, thickness=4], a=1..%, frames=180, background=P, size=[700,500], paraminfo=false);

                                      

                   

In the final example, we search for Hamiltonian paths in a lattice defined on the surface of a cube. Imagine a cube made of wires, and an ant must crawl along these wires from point  A(0,0,0)  to point B(2,2,2) , visiting all nodes of this lattice. Is this possible? We see that it is not. The length of the maximum path is 25, and this lattice has 26 vertices. An animation of one of the maximum paths is provided.

                           

S:={seq(seq(seq([i,j,k],i=0..2),j=0..2),k=0..2)} minus {[1,1,1]}: A:=[0,0,0]: B:=[2,2,2]:
P:=plots:-display(plots:-pointplot3d(S, symbol=solidcircle, color=blue, symbolsize=20, scaling=constrained), plots:-textplot3d([[A[],"A"],[B[],"B"]], font=[times,bold,22], align=[left,above]), plottools:-curve([[2,0,1],[2,2,1],[0,2,1],[0,0,1],[2,0,1]], color=black,thickness=0),plottools:-curve([[1,0,2],[1,0,0],[1,2,0],[1,2,2],[1,0,2]], color=black, thickness=0),plottools:-curve([[2,1,0],[0,1,0],[0,1,2],[2,1,2],[2,1,0]], color=black,thickness=0), tickmarks=[3,3,3]):
L:=AllPaths(S,A,B):
nops(L);
map(t->nops(t), L);
L1:=select(t->nops(t)=%[-1], L):
nops(L1);
plots:-animate((plots:-display)@(plottools:-curve),[L1[-1][1..round(a)], color=red, thickness=4], a=1..25, frames=240, background=P, paraminfo=false, axes=box, labels=[x,y,z]);

       

                              

We can see from this animation that the path does not pass through the vertex (0, 2, 2) .

Edit. A code error that could cause incorrect operation when using the  R  option has been fixed. Everything now works correctly. If there is no Hamiltonian path passing through all vertices, the procedure returns the empty set { } .

Paths1.mw

This system of nonlinear equations was proposed by an artificial intelligence named Alice.

f1 := x4^3+x5^2*x6-sin(x3)+exp(x1*x2)-5; 
f2 := ln(x4+x5)-x6^2+x3^3-x1^2*x2-2; 
f3 := cos(x4*x5)+x6*sin(x3)-x1*x2^2+x4^2*x5^3;


At the same time, AI accompanied it with comments. Here are some of them:
 

### System Characteristics
 **Underdetermination**: 6 unknowns with 3 equations 

### Solution Complexity

* Lack of a general analytical solution
* Need to use numerical methods
* Possibility of multiple solutions
* Complexity of visualization
* Sensitivity to initial approximations

### Solution Methods

* **Iterative methods**:
* Newton's method
* Simple iteration method
* Gradient methods
* **Numerical methods**:
* Finite difference method
* Monte Carlo method
* Genetic algorithms
* **Optimization approaches**:
* Lagrange multiplier method
* Constrained optimization methods

### Practical Application

Similar systems are found in:
* Quantum mechanics
* Field theory
* Economic modeling
* Bioinformatics
* Machine learning
* Engineering calculations of complex systems

Analyzing such systems often requires specialized software and powerful computing resources.

 


Although the worksheets provided here have been developed under Maple 2015, they should work correctly with newer versions, except perhaps for commands that use the 'op' function ('piecewise' mainly).

In the sequel the acronym 'pdf' stands for 'probability density function'.



CONTEXT

This post originates from a recent question by @JoyDivisionMan and the ensuing discussion. 

In a few words, the OP noticed Maple 2025 failed to return a result and asked why. In his reply, @acer identified a code regression somewhere in between Maple 2023 and Maple 2025.
Like Maple 2023 "my"  Maple 2015 does not fail but provide... a wrong answer: Do versions in between 2015 and 2024 return wrong results too?

 

In this post I explain how we can calculate the result by hand (only elementary maths required), why Maple 2015 (and likely newer versions) returns an incorrect result, why Maple generally fails in returning a result, and finally provide several examples to illustrate that even mathematically simple.
The various test cases I present are all equally simple for a skilled human agent, but conversely all beyond the reach of Statistics:-PDF.
This raises the question: Can a robust algorithm for calculating this PDF be developed without resorting to an expert system (I don't like the term AI) that mimics human reasoning?



THE MAIN OBSERVATION

Let X some continuous univariate random variable (CURV) and 𝟇 a real valued function from defined over the support of  . Let Y the random variable defined by Y =  𝟇(X).

The main observation about Statistics:-PDF is that unless very specific situations, this procedure does not build the correct expression of pdf
(Y) as soon as 𝟇 is not a strictly monotone function

Here are a fex exceptions to this claim:

X any CURV𝟇 : xxn  , n positive integer (correct solution even if 𝟇 is not monotone)

X ~ Uniform(-1, 1)𝟇 : xarctanh(x)  (no result returned even if 𝟇 is strictly monotone).


It is worth saying that, only by chance, Statistics:-PDF may sometimes return the correct result even if 𝟇 is a non monotone function.
For instance, in
@JoyDivisionMan's question, the procedure was indeed capable to provide a correct result for the case 
X ~ Uniform(0, 2𝜋)𝟇 : x ⟼ arctanh(x) but this was only because two errors balanced each other out. Replace X ~ Uniform(0, 2𝜋)  by X ~ Uniform(a, a+2𝜋) and Maple is wrong (notice that the pdf of 𝟇(X) remains unchanged whatever the value of a).


A GOOD DRAWING WORTH A THOUSAND WORDS

Here is a picture to help understand how to get the pdf of Y =  𝟇(X) for a non monotone function 𝟇 (the case of a monotone function directly comes from this latter).

In this illustration 
X ~ Uniform(0, 2𝜋)𝟇 : xsine(x).
To ease the explanation, I write 
X as a mixture of three uniform random variables X1, X2, X3, whose supports are the intervals of the three branches of 𝟇. More formally, X  = (1/4)∙X1 + (1/2)∙X2 + (1/4)∙X3.
The restrictions  of
𝟇 to these three branches are denoted 𝟇1, 𝟇2, 𝟇3.



The large rectangles below the horizontal axis represent the pdf of X1, X2, X3  and the blue curve the 𝟇 function.
The image of the 
interval [y-dy, y+dy]  by the inverse functions 𝟇1(-1) and 𝟇2(-1) of 𝟇1 and 𝟇2  are represented by the vertical rectangles [x1-dy, x1+dy] and [x2-dy, x2+dy] .
These two intervals bring a contribution to the pdf of 
(light gray blue on the right) represented by the horizontal violet rectangle on the right side of the picture.

T
he probability Prob(Y  [y-dy, y+dy]) that Y belongs to the interval [y-dy, y+dy]) is simply the sum 
                 Prob(
Y  [y-dy, y+dy])  = Prob(X1  𝟇1(-1)([y-dy, y+dy]))  + Prob(X2  𝟇2(-1)([y-dy, y+dy]))

Let
𝟇'b(-1)(y) denote the derivative of 𝟇b(-1)(y).
Making 
dy tends to 0 gives
                pdf(
Y y)  =  pdf(X1 𝟇1(-1)(y)) ×  𝟇'1(-1)(y) |pdf(X2 𝟇2(-1)(y)) ×  𝟇'2(-1)(y) |

As I said before there is truly no big math behind this, Except maybe those absolute values?
To understand where they come from zoom in on the rectangle 𝜔 = [
x1-dx, x1+dx] ╳ [y-dy, y+dy] and denote X𝜔 and Y𝜔 the restrictions of X1 and Y to 𝜔.
Locally 
Y𝜔 is proportional to A+B∙X𝜔 where constant B = 𝟇'(x1) and the value of constant does not matter here.
So the pdf of 
Y𝜔 is (a classical result)  pdf(Y𝜔 y) = pdf(X𝜔 (y-A)/C) / |C|.
Few details can be found Here.


MAPLE FAILURES AND WEAKNESSES

So why did Maple, at last some versions, produce a wrong result and why some versions are not even capable to return one?
The reason is that there is no big math only at first sight...because determining the inverse function of 
𝟇 can be quite tricky as soon as 𝟇 is not one-to-one map, for instance when 𝟇 is not strictly monotone.
When it is so
 𝟇(-1) must be defined for all the branches whose definition intervals intersect the support of X.

I spent a lot of time debugging the procedure Statistics:-PDF to understand why it either fails or produces incorrec results.
The 
sine_debug_nodebugoutput.mw  worksheet presents the "X ~ Uniform(0, 2𝜋)𝟇 : xsine(x)" case (as Mapleprimes stubbornly refuses to upload the worksheet containing the debugger trace, I convert it to sine_debug.pdf to help you see this trace). 
To orient the core development team correcting this procedure (assumming they care), the critical procedures are

Statistics:-RandomVariables:-PDF:-Univariate:-GetValueTab[anything]
and  Statistics:-RandomVariables:-GetInverse
 

At the very end it is this second procedure which is truly responsible of Maple failing to provide a result or returning an incorrect on, because it does not correctly build the inverse functions 𝟇b(-1)  for all the branches b which matter.
 

I wrote above that "it is only by chance that Maple provided the correct result for the non monotone function 𝟇 = cosine". Indeed Statistics:-PDF returns a wrong result when X ~ Uniform(z, z+2𝜋) and z is not a multiple of 𝜋/2 (see cosine.mw).

Other important situations where Maple fails returning a result are those where
𝟇 is a polynomial function with different zeros located in the support of X.
I did not trace them but it seems that
Statistics:-PDF does not know how to build the 𝟇b(-1) in this case (even though it is quite simple, see "Polynomial" examples below).


A SELECTION OF EXAMPLES

Here is a selection of examples to demonstrate that even in rather complex cases the pdfs of 𝟇(X) can be constructed quite easily (note that Maple either fails to compute them or to provide a correct result):



OPEN QUESTION

Tracing Statistics:-PDF reveals an already complex algorithm designed to handle a broad variety of "canonical" situations. Even this the algorithm fails in almost all non-toy-problem such as this compilation proves Maple_failures.mw.
At first sight, there seems to be a contradiction between the (apparent?) simplicity with which one can obtain, by hand in sometimes in an ad hoc way, the expression of pdf(𝟇(
X)), whether it is exact or truncated, and the complexity of the Statistics:-PDF algorithm, which results in failure in all non trivial cases.

This observation leads to the important question "Is it possible to rewrite Statistics:-PDF in order to enlarge its domain of success?".
I have the feeling that this means designing an algorithm which focuses more on mimicing the human reasoning than identifying "canonical" situations (as it is done today). 
An AI-driven algorithm maybe?

I repeat here that the maths are very simple, and all the more simple if you represent the random variable 
X as a mixture of components X1, ... XBboth having the same (truncated) distribution than X, and whose supports identify to the intervals of definition of the B branches of 𝟇 over the whole support of X.
The only difficulty lies in the identification of these branches and in the construction of the functions 
𝟇b(-1) over the supports of each Xb.

I have no answer to this question.

Hi there, distinguished audience,

Let f(n)=n2+n+41.

further, require n to be a positive integer.
Look at the cases when f(n) is a prime number, or a composite number.
Beautiful parabola patterns appear.
see attached

writeup_for_prime_producing_trinomial_final.pdf

I also look at the case when there is

n2+n+17.

see

Lucky Number of Euler -- from Wolfram MathWorld

I can be reached by email at

matthewcharlesanderson2@gmail.com

also,

https://mattanderson.fun

Regards,

Matthew

This post stems from this Question to which the author has never taken the time to give any answer whatsoever.

To help the reader understand what this is all about, I reproduce an abriged version of this question

I have the following data ... [and I want to]  create a cumulative histogram with corresponding polygon employing this same information...

The data the author refers to is a collection of decimal numbers.

The term "histogram" has a very well meaning in Statistics, without entering into technical details, let us say an histogram is an estimator of a Probability Density Function (continuous random variable) or of a mass function (discrete random variable), see for instance Freedman & Diaconis.

The expression "cumulative histogram" is more recent, see for instance Wiki for a quick explanation. Shortly a cumulative histogram can be seen as an approximation of the Cumulative Density Function (CDF) of the random variable whose the sample at hand is drawn from.

In fact there exists an alternative concept named ECDF (Empirical Cumulative Distribution Function) which has been around for a long time and which is already an estimator of the CDF.
Personally I am always surprised, given the many parameters it depends upon (anchors, number of bins, binwidth selection method, ...), when someone wants to draw a cumulative histogram: Why not draw instead the ECDF, a more objective estimator, even simpler to build than the cumulative histogram, and which does not use any parameter (that people often tune to get a pretty image instead of having a reliable estimator)? 

Anyway, I have done a little bit of work arround the OP's question, and it ended in a procedure named Hodgepodge (surely not a very explicit name but I was lacking inspiration) which enables plotting (if asked) several informations in addition to the required cumulative histogram:

  • The histogram of the raw data for the same list of bin bounds.
  • The kernel density estimator of this raw-data-histogram.
  • The ECDF of the data.

Here is an example of data

and here is what procedure Hodgepodge.mw  can display when all the graphics are requested

The inscribed square problem, also known as the Toeplitz conjecture, is an unsolved quastion in geometry: Does every plane simple closed curve (Jordan curve) contain all four vertices of some square? This is true if the curve is convex or piecewise smooth and in other special cases. The problem was proposed by Otto Toeplitz in 1911. For detailes see  https://en.wikipedia.org/wiki/Inscribed_square_problem

The Inscribed_Square procedure finds numerically one or more solutions for a curve defined by parametric equations of its boundary or by the equation F(x,y)=0. The required parameter of procedure  L  is the list of equations of the boundary links or the equation  F(x,y)=0 . Optional parameters:  N  and  R . By default  N='onesolution' (the procedure finds one solution), if  N  is any symbol (for example  N='s'), then more solutions.  R  is the range for the length of the side of the square (by defalt  R=0.1..100 ).

The second procedure  Pic  visualizes the results obtained.

The codes of the procedures:

restart;
Inscribed_Square:=proc(L::{list(list),`=`},N::symbol:='onesolution',R::range:=0.1..100)
local D, n, c, L1, L2, L3, f, L0, i, j, k, m, A, B, C, P, M, eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, sol, Sol;
uses LinearAlgebra;
if L::list then
L0:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L);
c:=0;
n:=nops(L);
for i from 1 to n do
for j from i to n do
for k from j to n do
for m from k to n do
A:=convert(subs(t=t1,L0[i,1]),Vector): 
B:=convert(subs(t=t2,L0[j,1]),Vector):
C:=convert(subs(t=t3,L0[k,1]),Vector): 
D:=convert(subs(t=t4,L0[m,1]),Vector):
M:=<0,-1;1,0>;
eq1:=eval(C[1])=eval((B+M.(B-A))[1]);
eq2:=eval(C[2])=eval((B+M.(B-A))[2]);
eq3:=eval(D[1])=eval((C+M.(C-B))[1]);
eq4:=eval(D[2])=eval((C+M.(C-B))[2]);
eq5:=eval(DotProduct(B-A,B-A, conjugate=false))=d^2;
sol:=fsolve([eq1,eq2,eq3,eq4,eq5],{t1=op([2,2,1],L0[i])..op([2,2,2],L0[i]),t2=op([2,2,1],L0[j])..op([2,2,2],L0[j]),t3=op([2,2,1],L0[k])..op([2,2,2],L0[k]),t4=op([2,2,1],L0[m])..op([2,2,2],L0[m]),d=R});
if type(sol,set(`=`)) then if N='onesolution' then return convert~(eval([A,B,C,D],sol),list) else c:=c+1; Sol[c]:=convert~(eval([A,B,C,D],sol),list) fi;
 fi; 
od: od: od: od:
Sol:=fnormal(convert(Sol,list),7);
print(Sol);
ListTools:-Categorize((X,Y)->`and`(seq(is(convert(X,set)[i]=convert(Y,set)[i]),i=1..4)) , Sol);
return map(t->t[1],[%]);
else
A,B,C,D:=<x1,y1>,<x2,y2>,<x3,y3>,<x4,y4>:
M:=<0,-1;1,0>:
eq1:=eval(C[1])=eval((B+M.(B-A))[1]):
eq2:=eval(C[2])=eval((B+M.(B-A))[2]):
eq3:=eval(D[1])=eval((C+M.(C-B))[1]):
eq4:=eval(D[2])=eval((C+M.(C-B))[2]):
eq5:=eval(LinearAlgebra:-DotProduct((B-A,B-A), conjugate=false))=d^2:
eq6:=eval(L,[x=x1,y=y1]):
eq7:=eval(L,[x=x2,y=y2]):
eq8:=eval(L,[x=x3,y=y3]):
eq9:=eval(L,[x=x4,y=y4]):
sol:=fsolve({eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9},{seq([x||i=-2..2,y||i=-2..2][],i=1..4),d=R});
eval([[x1,y1],[x2,y2],[x3,y3],[x4,y4]], sol):
fi;
end proc:

Pic:=proc(L,Sol,R::range:=-20..20)
local P1, P2, P3, T;
uses plots, plottools;
P1:=`if`(L::list,seq(`if`(type(s,listlist),line(s[],color=blue, thickness=2),plot([s[1][],s[2]],color=blue, thickness=2)),s=L), implicitplot(L, x=R,y=R, color=blue, thickness=2, gridrefine=3));
P2:=polygon(Sol,color=yellow,thickness=0);
P3:=curve([Sol[],Sol[1]],color=red,thickness=3):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"]], font=[times,18], align=[left,above]);
display(P1,P2,P3,T, scaling=constrained, size=[800,500], axes=none);
end proc:

Examples of use:

The curve consists of a semicircle, a segment and a semi-ellipse (find 1 solution):

L:=[[[cos(t),sin(t)],t=0..Pi],[[t,0],t=-1..0],[[0.5+0.5*cos(t),0.8*sin(t)],t=Pi..2*Pi]]:
Sol:=Inscribed_Square(L);
Pic(L,Sol);

       


The procedure finds 6 solutions for a non-convex pentagon:

 L:=[[[0,0],[9,0]],[[9,0],[8,5]],[[8,5],[5,3]],[[5,3],[0,4]],[[0,4],[0,0]]]:
Sol:=Inscribed_Square(L,'s');
plots:-display(Matrix(3,2,[seq(Pic(L,Sol[i]),i=1..6)]),size=[300,200]);

             


For an implicitly defined curve, only one solution can be found:

L:=abs(x)+2*abs(y)-sin((2*x-y))-cos(x+y)^2=3:
Sol:=Inscribed_Square(L);
Pic(L,Sol);

               
See more examples in the attached file.

Inscribed_Square.mw

 

Hi again all,

prime numbers are fun for me.

see

pairs_of_prime_numbers_procedure_with_union_and_isprime.pdf

sorry, could not find the .mw file,

but the code is small, and easy to copy

this is fun for me, here in Keizer OR, USA

Party on, everyone

Best regards,

Matt

https://mattanderson.fun/

My student, David Wei, has ported Andrew Hicks's MOISE package to Maple 2025:

GitHub - david-wei-01001/MOISE-for-Maple-2025: MOISE that is compatible with the newest Maple 2025 Engine

Regards,

Gerald Penn

The Autumn Issue is now up, at mapletransactions.org

This issue contains two Featured Contributions; a short but very interesting one by Gilbert Labelle on a topic very dear to my own heart, and a longer and also very interesting one by Wadim Zudilin.  I asked Doron Zeilberger about Wadim's paper, and he said "this is a true gem with lots of insight and making connections between different approaches."

The "Editor's Corner" paper is a little different, this time.  This paper is largely the work of my co-author, Michelle Hatzel, extracted and revised from her Masters' thesis which she defended successfully this past August.  I hope that you find it as interesting as I did.

 

We have three refereed contributions, a contribution on the use of Maple Learn in teaching, and a little note on my design of the 2026 Calendar for my upcoming SIAM book with Nic Fillion, as well.  All the images for the calendar were generated in Maple (as were most of the images in the book).

It's been fun to put this issue together (with an enormous amount of help from Michelle) and I hope that you enjoy reading it.

I would also like to thank the Associate Editors who handled the refereeing: Dhavide Aruliah, David Jeffrey, and Viktor Levandovskyy.

It is possible to construct a system of equations that defines all spheres that touch two smooth surfaces.

This includes all the spheres inscribed between these surfaces. The first two equations are surfaces,  x7, x8, x9 are the coordinates of the centers of the spheres, the remaining variables correspond to the coordinates of the points of contact. 
In the program, we inscribe a sphere between a cylinder and a sphere. An example of obtaining one of the infinite subsets of the solution to a system of equations is implemented. This solution is "responsible" for the equation f8.
INSCRIBED_SPHERES_CILL_FOR.mw

And one more example.

[Edit: a general procedure based on these ideas is given below.]

Maple's isolve can solve some but not all of the quadratic equations in two variables describing the conic sections. Here I show that by transforming the equations, Maple can then solve these missing cases.

restart

with(plots)

I was recently surprised that isolve cannot solve the following simple Diophantine equation

isolve(2*x^2+4*x*y+y^2 = 7)

which has the obvious (to a human) solution {x = 1, y = 1}. This led me to think about the case of conic sections, which have the following general equation (= 0 implied), where I assume a, () .. (), f are integers.

P := a*x^2+b*x*y+c*y^2+d*x+e*y+f

a*x^2+b*x*y+c*y^2+d*x+e*y+f

The above equation has discriminant -4*a*c+b^2 positive, indicating that is is a hyperbola.

h_params := {a = 2, b = 4, c = 1, d = 0, e = 0, f = -7}; P__h := eval(P, h_params); disc := -4*a*c+b^2; eval(disc, h_params); plot__h := implicitplot(P__h, x = -10 .. 10, y = -10 .. 10, colour = red)

{a = 2, b = 4, c = 1, d = 0, e = 0, f = -7}

2*x^2+4*x*y+y^2-7

-4*a*c+b^2

8

Here's a parabola case (discriminant = 0) that isolve also has trouble with

p_params := {a = 2, b = 4, c = 2, d = 1, e = 2, f = 7}; eval(disc, p_params); P__p := eval(P, p_params); {isolve(P__p)}; plot__h := implicitplot(P__p, x = -10 .. 10, y = -10 .. 10, colour = red)

{a = 2, b = 4, c = 2, d = 1, e = 2, f = 7}

0

2*x^2+4*x*y+2*y^2+x+2*y+7

{}

But this has at least one solution

eval(P__p, {x = 7, y = -7})

0

Maple seems to do better in the elliptic case (discriminant negative) and finds two solutions. Examination of the plot suggests there are no other solutions.

e_params := {a = 2, b = 4, c = 3, d = 1, e = 2, f = -7}; eval(disc, e_params); P__e := eval(P, e_params); {isolve(P__e)}; plot__e := implicitplot(P__e, x = -10 .. 10, y = -10 .. 10, colour = red)

{a = 2, b = 4, c = 3, d = 1, e = 2, f = -7}

-8

2*x^2+4*x*y+3*y^2+x+2*y-7

{{x = -3, y = 2}, {x = 2, y = -3}}

I show that transformation of the general equation to the form -D*Y^2+X^2 = m, where D is the discriminant, allows Maple to solve the hyperbolic case, as well as the elliptic case it already knows how to solve; another transformation works for the parabolic case. Maple appears to be able to solve all (solvable) cases of the transformed equations, though this is not clear from the help page. The transformation is discussed in Bogdan Grechuk, Polynomial Diophantine Equations: A Systematic Approach, Springer, 2024, Sec. 3.1.7. doi: 10.1007/978-3-031-62949-5

 

A complete classification of the conics, including degenerate cases, is given in https://en.wikipedia.org/wiki/Matrix_representation_of_conic_sections. If the determinant, delta, of the following matrix is zero, we have a degenerate case.

A__Q := Matrix(3, 3, [a, (1/2)*b, (1/2)*d, (1/2)*b, c, (1/2)*e, (1/2)*d, (1/2)*e, f]); delta := LinearAlgebra:-Determinant(A__Q)

Matrix(%id = 36893489963432522084)

a*c*f-(1/4)*a*e^2-(1/4)*b^2*f+(1/4)*b*d*e-(1/4)*c*d^2

The case of a = b and b = c and c = 0 is just the linear case, which Maple can solve, and is one case where D = delta and delta = 0. Other degenerate parabola cases are two coincident lines or two parallel lines. Degenerate hyperbolas are two intersecting lines, and degenerate ellipses are a single point. Maple can solve all these cases, e.g.,NULL

expand((2*x+3*y+1)*(2*x+3*y)); {isolve(%)}; expand((2*x+3*y)*(2*x+3*y)); {isolve(%)}; expand((x-2)*(y-3)); {isolve(%)}; x^2+y^2 = 0; {isolve(%)}

4*x^2+12*x*y+9*y^2+2*x+3*y

{{x = -3*_Z1, y = 2*_Z1}, {x = -2-3*_Z1, y = 1+2*_Z1}}

4*x^2+12*x*y+9*y^2

{{x = -3*_Z1, y = 2*_Z1}}

x*y-3*x-2*y+6

{{x = _Z1, y = 3}}

x^2+y^2 = 0

{{x = 0, y = 0}}

(The intersecting lines case above only finds one of the lines.) The transformation will consider only the case where at least one of a or c is non-zero. This misses hyperbolas with a = c and c = 0; Maple seems to handle these bilinear equations, e.g.,

bl_params := {a = 0, b = 2, c = 0, d = 1, e = 1, f = 2}; P__bl := eval(P, bl_params); {isolve(P__bl)}

{a = 0, b = 2, c = 0, d = 1, e = 1, f = 2}

2*x*y+x+y+2

{{x = -2, y = 0}, {x = -1, y = 1}, {x = 0, y = -2}, {x = 1, y = -1}}

Transformation for non-zero discriminant
At least one of a or c must be non-zero; if necessary exchange x and y to ensure a is non-zero.

We multiply P by -4*a*(-4*a*c+b^2) and change to new variables X and Y.

itr := {X = (-4*a*c+b^2)*y+b*d-2*a*e, Y = 2*a*x+b*y+d}; tr := solve(itr, {x, y})

{X = (-4*a*c+b^2)*y+b*d-2*a*e, Y = 2*a*x+b*y+d}

{x = (1/2)*(4*Y*a*c-Y*b^2+2*a*b*e-4*a*c*d+X*b)/((4*a*c-b^2)*a), y = -(2*a*e-b*d+X)/(4*a*c-b^2)}

The transformed equation has the form -D*Y^2+X^2-m = 0 or -D*Y^2+X^2 = m, where D is the discriminant.

Q := collect(normal(eval(-4*P*a*(-4*a*c+b^2), tr)), {X, Y}); m := -coeff(coeff(Q, X, 0), Y, 0)

X^2+(4*a*c-b^2)*Y^2+16*a^2*c*f-4*a^2*e^2-4*a*b^2*f+4*a*b*d*e-4*a*c*d^2

-16*a^2*c*f+4*a^2*e^2+4*a*b^2*f-4*a*b*d*e+4*a*c*d^2

For positive discriminant D, this is a general Pell's equation, -D*Y^2+X^2 = m, which Maple knows how to solve. (The definition of Pell's equation requires that D is not a square, but Maple can also solve the simpler case where D is a square.) For the hyperbola above, we have an infinite number of solutions, parameterized by an arbitrary integer _Z1.

Q__h := eval(Q, h_params); solXY__h := {isolve(Q__h)}

X^2-8*Y^2+448

{{X = -12*(1+2^(1/2))^(1+2*_Z1)-12*(1-2^(1/2))^(1+2*_Z1)-4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = -3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))-2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)}, {X = -12*(1+2^(1/2))^(1+2*_Z1)-12*(1-2^(1/2))^(1+2*_Z1)-4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = 3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))+2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)}, {X = -12*(1+2^(1/2))^(1+2*_Z1)-12*(1-2^(1/2))^(1+2*_Z1)+4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = -3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))+2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)}, {X = -12*(1+2^(1/2))^(1+2*_Z1)-12*(1-2^(1/2))^(1+2*_Z1)+4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = 3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))-2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)}, {X = 12*(1+2^(1/2))^(1+2*_Z1)+12*(1-2^(1/2))^(1+2*_Z1)-4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = -3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))+2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)}, {X = 12*(1+2^(1/2))^(1+2*_Z1)+12*(1-2^(1/2))^(1+2*_Z1)-4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = 3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))-2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)}, {X = 12*(1+2^(1/2))^(1+2*_Z1)+12*(1-2^(1/2))^(1+2*_Z1)+4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = -3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))-2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)}, {X = 12*(1+2^(1/2))^(1+2*_Z1)+12*(1-2^(1/2))^(1+2*_Z1)+4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = 3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))+2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)}}

Transforming back to the original coordinates

sol__h := {seq(eval(eval(tr, h_params), solXY), `in`(solXY, solXY__h))}

{{x = -2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)-(5/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = (3/2)*(1+2^(1/2))^(1+2*_Z1)+(3/2)*(1-2^(1/2))^(1+2*_Z1)+(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = -2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)+(5/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = (3/2)*(1+2^(1/2))^(1+2*_Z1)+(3/2)*(1-2^(1/2))^(1+2*_Z1)-(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = -(1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)-(1/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = (3/2)*(1+2^(1/2))^(1+2*_Z1)+(3/2)*(1-2^(1/2))^(1+2*_Z1)-(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = -(1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)+(1/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = (3/2)*(1+2^(1/2))^(1+2*_Z1)+(3/2)*(1-2^(1/2))^(1+2*_Z1)+(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = (1+2^(1/2))^(1+2*_Z1)+(1-2^(1/2))^(1+2*_Z1)-(1/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = -(3/2)*(1+2^(1/2))^(1+2*_Z1)-(3/2)*(1-2^(1/2))^(1+2*_Z1)-(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = (1+2^(1/2))^(1+2*_Z1)+(1-2^(1/2))^(1+2*_Z1)+(1/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = -(3/2)*(1+2^(1/2))^(1+2*_Z1)-(3/2)*(1-2^(1/2))^(1+2*_Z1)+(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = 2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)-(5/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = -(3/2)*(1+2^(1/2))^(1+2*_Z1)-(3/2)*(1-2^(1/2))^(1+2*_Z1)+(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = 2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)+(5/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = -(3/2)*(1+2^(1/2))^(1+2*_Z1)-(3/2)*(1-2^(1/2))^(1+2*_Z1)-(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}}

Evaluate for some _Z1 values, say _Z1=0 and _Z1=5.
It is evident from tr above that integer solutions in X, Y may transform to non-integer solutions in x, y but that doesn't occur for these two cases

sol__h0 := evala(eval(sol__h, _Z1 = 0)); sol__h5 := evala(eval(sol__h, _Z1 = 5))

{{x = -9, y = 5}, {x = -3, y = 1}, {x = -1, y = -1}, {x = -1, y = 5}, {x = 1, y = -5}, {x = 1, y = 1}, {x = 3, y = -1}, {x = 9, y = -5}}

{{x = -61181, y = 35839}, {x = -21979, y = 12875}, {x = -10497, y = 35839}, {x = -3771, y = 12875}, {x = 3771, y = -12875}, {x = 10497, y = -35839}, {x = 21979, y = -12875}, {x = 61181, y = -35839}}

Check they are solutions to P__h

map2(eval, P__h, `union`(sol__h0, sol__h5))

{0}

For negative discriminant and negative m, the equation -D*Y^2+X^2 = m has no solutions. In the classification scheme for the conics, the "imaginary ellipse" case (no real solutions) occurs when (a+c)*delta > 0. For negative discriminant, we must have a and c the same sign, and this is the case of negative m.

factor(expand((16*(a+c))*delta)); factor(-m)

4*(a+c)*(4*a*c*f-a*e^2-b^2*f+b*d*e-c*d^2)

4*a*(4*a*c*f-a*e^2-b^2*f+b*d*e-c*d^2)

. In this case Maple returns NULL for both the untransformed and transformed case, which can mean no solutions or just that isolve couldn't find any.

ie_params := {a = 3, b = 0, c = 5, d = 0, e = 0, f = 8}; P__ie := eval(P, ie_params); {isolve(P__ie)}; Q__ie := eval(Q, ie_params); {isolve(Q__ie)}

{a = 3, b = 0, c = 5, d = 0, e = 0, f = 8}

3*x^2+5*y^2+8

{}

X^2+60*Y^2+5760

{}

For negative discriminant and positive m or (a+c)*delta < 0, the ellipse is real, and there are a finite number of solutions. Maple solves the untransformed and transformed equations.
Here we need to filter out non-integer solutions

P__e; {isolve(P__e)}; Q__e := eval(Q, e_params); solXY__e := {isolve(Q__e)}; {seq(eval(eval(tr, e_params), solXY), `in`(solXY, solXY__e))}; sol__e := select(proc (q) options operator, arrow; eval(x::integer, q) and eval(y::integer, q) end proc, %)

2*x^2+4*x*y+3*y^2+x+2*y-7

{{x = -3, y = 2}, {x = 2, y = -3}}

X^2+8*Y^2-472

{{X = -20, Y = -3}, {X = -20, Y = 3}, {X = 20, Y = -3}, {X = 20, Y = 3}}

{{x = -3, y = 2}, {x = 2, y = -3}, {x = -3/2, y = 2}, {x = 7/2, y = -3}}

{{x = -3, y = 2}, {x = 2, y = -3}}

Transformation for zero discriminant
For the case of zero discriminant (parabola), we need a different transformation.

itr0 := {X = 2*a*x+b*y+d, Y = y}; tr0 := solve(itr0, {x, y})

{X = 2*a*x+b*y+d, Y = y}

{x = (1/2)*(-Y*b+X-d)/a, y = Y}

The transformed equation is of the form A*Y+X^2+B = 0

Q0 := collect(normal(eval(4*a*P, tr0)), {X, Y})

X^2+(4*a*c-b^2)*Y^2+(4*a*e-2*b*d)*Y+4*f*a-d^2

We consider the parabolic example above, for which Maple finds no solutions without transformation.

P__p; {isolve(P__p)}

2*x^2+4*x*y+2*y^2+x+2*y+7

{}

For the transformed problem, Maple finds an infinite number of solutions

Q__p := eval(Q0, p_params); solXY__p := {isolve(Q__p)}; sol__p := {seq(eval(eval(tr0, p_params), solXY), `in`(solXY, solXY__p))}

X^2+8*Y+55

{{X = 1+8*_Z1, Y = -8*_Z1^2-2*_Z1-7}, {X = 3+8*_Z1, Y = -8*_Z1^2-6*_Z1-8}, {X = 5+8*_Z1, Y = -8*_Z1^2-10*_Z1-10}, {X = 7+8*_Z1, Y = -8*_Z1^2-14*_Z1-13}}

{{x = 17/2+8*_Z1+8*_Z1^2, y = -8*_Z1^2-6*_Z1-8}, {x = 29/2+16*_Z1+8*_Z1^2, y = -8*_Z1^2-14*_Z1-13}, {x = 8*_Z1^2+4*_Z1+7, y = -8*_Z1^2-2*_Z1-7}, {x = 8*_Z1^2+12*_Z1+11, y = -8*_Z1^2-10*_Z1-10}}

Two of the general solutions will not give integer solutions, so could be filtered out, but it is easier to filter after choosing some specific _Z1 values.

eval(sol__p, _Z1 = 0); sol__p0 := select(proc (q) options operator, arrow; eval(x::integer, q) and eval(y::integer, q) end proc, %); eval(sol__p, _Z1 = 5); sol__p5 := select(proc (q) options operator, arrow; eval(x::integer, q) and eval(y::integer, q) end proc, %)

{{x = 7, y = -7}, {x = 11, y = -10}, {x = 17/2, y = -8}, {x = 29/2, y = -13}}

{{x = 7, y = -7}, {x = 11, y = -10}}

{{x = 227, y = -217}, {x = 271, y = -260}, {x = 497/2, y = -238}, {x = 589/2, y = -283}}

{{x = 227, y = -217}, {x = 271, y = -260}}

Check they are solutions to P__p

map2(eval, P__p, `union`(sol__p0, sol__p5))

{0}

NULL

Download DiophantineConics2.mw

It's possible to carve a hole through a unit cube, without splitting it into pieces, so that another unit cube can pass through that hole.  This is know as the Prince Rupert problem and was first analyzed by John Wallis, a contemporary of Isaac Newton.  Here's what the result looks like:

If your computer can play audio, have a look at  Ruperts Cube with music!

Here is the worksheet that produced the cube and the animation: ruperts-cube.mw

1 2 3 4 5 6 7 Last Page 1 of 78