#QUADRATIC DIOPHANTINE SOLVER 
# ==========================================================
# QUADRATIC DIOPHANTINE SOLVER - EDUCATIONAL ENGLISH EDITION
# Exact classification + integer search + global factor proof
# Load into Maple as a Maple input file (.mpl).
# ==========================================================
with(plots):
ClassifyQuadratic := proc(a::integer,b::integer,c::integer,
                          d::integer,e::integer,f::integer)
    local D, K4, label, center_value;

    D := b^2 - 4*a*c;
    # K4 is four times the determinant of the 3x3 conic matrix.
    # Exact integer arithmetic avoids floating-point comparisons.
    K4 := 4*a*c*f - a*e^2 - b^2*f + b*d*e - c*d^2;

    if a=0 and b=0 and c=0 then
        if d=0 and e=0 then
            if f=0 then
                label := "IDENTITY (ALL INTEGER PAIRS)";
            else
                label := "CONTRADICTION (NO SOLUTIONS)";
            end if;
        else
            label := "LINEAR EQUATION";
        end if;
    elif K4=0 then
        if D>0 then
            label := "DEGENERATE: TWO INTERSECTING LINES";
        elif D=0 then
            label := "DEGENERATE: PARALLEL/COINCIDENT OR NONREAL LINES";
        else
            label := "DEGENERATE: ONE REAL POINT";
        end if;
    elif D>0 then
        label := "HYPERBOLA (NON-DEGENERATE)";
    elif D=0 then
        label := "PARABOLA (NON-DEGENERATE)";
    else
        # For D<0 the quadratic form is definite.
        # The equation evaluated at the center is K4/(-D).
        center_value := K4/(-D);
        if a*center_value<0 then
            label := "REAL ELLIPSE (NON-DEGENERATE)";
        else
            label := "IMAGINARY ELLIPSE (NO REAL POINTS)";
        end if;
    end if;

    return table([
        "type"=label,
        "discriminant"=D,
        "determinant"=K4/4,
        "determinant_maal_4"=K4,
        "degenerate"=evalb(K4=0 and not (a=0 and b=0 and c=0))
    ]);
end proc:

FindIntegerSolutionsInBox := proc(a::integer,b::integer,c::integer,
                                  d::integer,e::integer,f::integer,
                                  xr::range,yr::range,limit::posint)
    local result, xx, yy, A, B, C, Dy, s, numer, denom;
    local xmin,xmax,ymin,ymax,sgn;

    xmin:=lhs(xr); xmax:=rhs(xr);
    ymin:=lhs(yr); ymax:=rhs(yr);
    result:={};

    for xx from xmin to xmax do
        A:=c;
        B:=b*xx+e;
        C:=a*xx^2+d*xx+f;

        if A<>0 then
            Dy:=B^2-4*A*C;
            if Dy>=0 and issqr(Dy) then
                s:=isqrt(Dy);
                denom:=2*A;
                for sgn in [-1,1] do
                    numer:=-B+sgn*s;
                    if irem(numer,denom)=0 then
                        yy:=iquo(numer,denom);
                        if yy>=ymin and yy<=ymax then
                            result:=result union {{x=xx,y=yy}};
                        end if;
                    end if;
                end do;
            end if;
        elif B<>0 then
            if irem(-C,B)=0 then
                yy:=iquo(-C,B);
                if yy>=ymin and yy<=ymax then
                    result:=result union {{x=xx,y=yy}};
                end if;
            end if;
        elif C=0 then
            # Iedere y in het zoekgebied is een oplossing.
            for yy from ymin to ymax do
                result:=result union {{x=xx,y=yy}};
                if nops(result)>=limit then
                    return result;
                end if;
            end do;
        end if;

        if nops(result)>=limit then
            return result;
        end if;
    end do;

    return result;
end proc:

PlotQuadraticFriendly := proc(P, solutions::set, xr::range, yr::range,
                              label::string)
    local curve, points, sol, xx, yy, layers, pt, pointLabel, plotTitle;
    local a0,b0,c0,d0,e0,den,xc,yc,vc,signQ,r,base,delta,j,L;
    local isImaginary, colors, centerPlot, firstLevel, targetLevels;

    points := [];
    for sol in solutions do
        if type(sol,set) then
            xx := eval(x,sol); yy := eval(y,sol);
            if type(xx,integer) and type(yy,integer)
               and xx>=lhs(xr) and xx<=rhs(xr)
               and yy>=lhs(yr) and yy<=rhs(yr) then
                points := [op(points),[xx,yy]];
            end if;
        end if;
    end do;
    points := sort(points,
        (u,v) -> evalb(u[1]<v[1] or (u[1]=v[1] and u[2]<v[2])));

    isImaginary := evalb(label="IMAGINARY ELLIPSE (NO REAL POINTS)");

    if isImaginary then
        # The original equation P=0 has no real curve.
        # Show nearby nonzero LEVEL SETS instead: these are real ellipses
        # and MUST NOT be described as solutions of P=0.
        a0 := coeff(P,x,2);
        b0 := coeff(coeff(P,x,1),y,1);
        c0 := coeff(P,y,2);
        d0 := eval(diff(P,x),{x=0,y=0});
        e0 := eval(diff(P,y),{x=0,y=0});
        den := 4*a0*c0-b0^2;
        xc := (b0*e0-2*c0*d0)/den;
        yc := (b0*d0-2*a0*e0)/den;
        vc := simplify(eval(P,{x=xc,y=yc}));
        signQ := signum(a0);
        # Begin at the next multiple of 10 beyond the extremum.
        # Example: P_min = 85/12 gives P = 10, 20, 40.
        firstLevel := 10*(floor(evalf(abs(vc)/10))+1);
        if firstLevel <= abs(vc) then
            firstLevel := firstLevel + 10;
        end if;
        targetLevels := [firstLevel,2*firstLevel,4*firstLevel];
        colors := [blue,green,magenta];
        layers := [];
        for j from 1 to 3 do
            L := signQ*targetLevels[j];
            try
                curve := plots:-implicitplot(P=L,x=xr,y=yr,
                         grid=[100,100],color=colors[j],thickness=2,
                         legend=sprintf("Level P = %a",L));
                layers := [op(layers),curve];
            catch:
            end try;
        end do;
        if xc>=lhs(xr) and xc<=rhs(xr)
           and yc>=lhs(yr) and yc<=rhs(yr) then
            if signQ=1 then
                pointLabel := sprintf("Minimum of P (%.2f, %.2f)",
                                      evalf(xc),evalf(yc));
            else
                pointLabel := sprintf("Maximum of P (%.2f, %.2f)",
                                      evalf(xc),evalf(yc));
            end if;
            centerPlot := plots:-pointplot([[xc,yc]],
                        symbol=solidcircle,symbolsize=14,color=red,
                        legend=pointLabel);
            layers := [op(layers),centerPlot];
        end if;
        plotTitle := "LEVEL CURVES: P=0 has no real solutions";
        return plots:-display(layers,axes=normal,
                    title=plotTitle,labels=["x","y"]);
    end if;

    try
        curve := plots:-implicitplot(P=0, x=xr, y=yr,
                    grid=[100,100], color=blue, thickness=2,
                    legend="Equation P(x,y) = 0");
    catch:
        curve := plot([],x=xr,y=yr);
    end try;

    layers := [curve];
    for pt in points do
        pointLabel := sprintf("(%a, %a)",pt[1],pt[2]);
        layers := [op(layers),
            plots:-pointplot([pt], color=red,
                symbol=solidcircle, symbolsize=14,
                legend=pointLabel)];
    end do;

    return plots:-display(layers, axes=normal, title=label,
                          labels=["x","y"]);
end proc:


# Exact finite divisor method for non-degenerate hyperbolas
# with square positive discriminant. Does not rely on a search box.
SquareHyperbolaProof := proc(a::integer,b::integer,c::integer,
                            d::integer,e::integer,f::integer)
    local A,B,C,E,F,G,D,k,h,N,pairs,u,v,dd,sgn,tval,yval,xval;
    local swapped,sols,divs;
    if a<>0 then
        A:=a; B:=b; C:=c; E:=d; F:=e; G:=f;
        swapped:=false;
    elif c<>0 then
        A:=c; B:=b; C:=a; E:=e; F:=d; G:=f;
        swapped:=true;
    else
        return table(["proved"=false,"solutions"={},
                      "reason"="Neither squared coefficient is nonzero."]);
    end if;
    D:=B^2-4*A*C;
    if D<=0 or not issqr(D) then
        return table(["proved"=false,"solutions"={},
                      "reason"="The positive discriminant is not a square."]);
    end if;
    k:=isqrt(D);
    h:=B*E-2*A*F;
    N:=D*(E^2-4*A*G)-h^2;
    if N=0 then
        return table(["proved"=false,"solutions"={},
                      "reason"="The conic is degenerate; a different proof is needed."]);
    end if;
    sols:={};
    divs:=numtheory:-divisors(abs(N));
    for dd in divs do
        for sgn in [-1,1] do
            u:=sgn*dd;
            v:=iquo(N,u);
            tval:=(u+v)/(2*k);
            yval:=((v-u)/2-h)/D;
            xval:=(tval-B*yval-E)/(2*A);
            if type(xval,integer) and type(yval,integer) then
                if swapped then
                    sols:=sols union {{x=yval,y=xval}};
                else
                    sols:=sols union {{x=xval,y=yval}};
                end if;
            end if;
        end do;
    end do;
    return table(["proved"=true,"solutions"=sols,
                  "D"=D,"k"=k,"h"=h,"N"=N,
                  "swapped"=swapped]);
end proc:

# Center and asymptotes exist for non-degenerate central conics.
CenterAndAsymptotes := proc(a::integer,b::integer,c::integer,
                            d::integer,e::integer)
    local D,xc,yc,X,Y,Q;
    D:=b^2-4*a*c;
    if D=0 then return NULL; end if;
    xc:=(2*c*d-b*e)/D;
    yc:=(2*a*e-b*d)/D;
    X:=x-xc; Y:=y-yc;
    Q:=expand(a*X^2+b*X*Y+c*Y^2);
    return table(["x"=xc,"y"=yc,"asymptotes"=factor(Q)]);
end proc:



# Friendly Level 2 explanation: algebra, finite factor pairs, and proof.
ExplainStudentLevel := proc(a::integer,b::integer,c::integer,
                            d::integer,e::integer,f::integer,
                            info, proof, results::set, complete::truefalse,
                            search_x::range, search_y::range)
    local D,K,A,B,C,E,F,G,k,h,N,T,U,V,sol,xx,yy,qq,divcount;
    D:=info["discriminant"];
    K:=info["determinant"];
    printf("\n=======================================================\n");
    printf("  LEARNING CORNER - LEVEL 2: STEP-BY-STEP EXPLANATION\n");
    printf("=======================================================\n");
    printf("STEP 1. Recognize the curve.\n");
    printf("  The discriminant is D = b^2 - 4*a*c = %a.\n",D);
    printf("  The conic determinant is K = %a.\n",K);
    if a=0 and b=0 and c=0 then
        printf("  There are no squared terms, so this is a linear/constant equation.\n");
    elif K=0 then
        printf("  K = 0: the conic is degenerate (not an ordinary smooth conic).\n");
    elif D>0 then
        printf("  D > 0 and K is not zero: the graph is a genuine hyperbola.\n");
    elif D=0 then
        printf("  D = 0 and K is not zero: the graph is a parabola.\n");
    else
        printf("  D < 0: the graph has elliptic type.\n");
        printf("  The determinant and coefficient signs decide whether real points exist.\n");
    end if;

    if complete then
        printf("\nSTEP 2. Replace the quadratic equation with a product.\n");
        if proof["swapped"] then
            printf("  For this calculation, X=y and Y=x (the variables are swapped).\n");
            A:=c; B:=b; C:=a; E:=e; F:=d; G:=f;
        else
            printf("  Here X=x and Y=y.\n");
            A:=a; B:=b; C:=c; E:=d; F:=e; G:=f;
        end if;
        k:=proof["k"]; h:=proof["h"]; N:=proof["N"];
        T:=2*A*'X'+B*'Y'+E;
        U:=expand(k*T-(D*'Y'+h));
        V:=expand(k*T+(D*'Y'+h));
        printf("  Define T = %a.\n",T);
        printf("  Define U = %a.\n",U);
        printf("  Define V = %a.\n",V);
        printf("  Every solution must satisfy U*V = %a.\n",N);
        if a=1 and b=8 and c=7 and d=2 and e=11 and f=12 then
            printf("  For this example we can divide both factors by 6:\n");
            printf("    (2*x + 2*y + 3)*(2*x + 14*y + 1) = -45\n");
            printf("  Put u=2*x+2*y+3 and v=2*x+14*y+1.\n");
            printf("  Then u*v=-45, which is easy to examine with divisors.\n");
        end if;
        printf("\nSTEP 3. Use the factor pairs.\n");
        printf("  An integer factor of %a must divide it exactly.\n",N);
        printf("  So there are only finitely many signed possibilities.\n");
        printf("  Test every divisor U of %a and set V = %a/U.\n",N,N);
        printf("  Number of signed possibilities to check: %d.\n",
               2*nops(numtheory:-divisors(abs(N))));
        printf("\nSTEP 4. Recover x and y, and keep integers only.\n");
        printf("  From U and V we get:\n");
        printf("    T = (U + V)/(2*%a)\n",k);
        printf("    Y = ((V - U)/2 - (%a))/%a\n",h,D);
        printf("    X = (T - (%a)*Y - (%a))/(2*(%a))\n",B,E,A);
        if proof["swapped"] then
            printf("  Remember that X=y and Y=x here.\n");
        end if;
        printf("  Fractions are rejected: both original coordinates must be integers.\n");
        printf("\nSTEP 5. State the answer and why it is complete.\n");
        printf("  Exactly %d distinct integer pairs satisfy the equation.\n",nops(results));
        printf("  All signed factor pairs were checked, not only points on the graph.\n");
        printf("  Therefore NO other integer solutions exist anywhere.\n");
    else
        printf("\nSTEP 2. Search for integer solutions.\n");
        printf("  For each integer x in %a, treat the equation as a quadratic in y.\n",search_x);
        printf("  A valid y must be an integer and lie in %a.\n",search_y);
        printf("\nSTEP 3. Interpret the result carefully.\n");
        printf("  %d distinct solutions were returned.\n",nops(results));
        printf("  This does NOT prove there are no solutions outside the search box.\n");
        printf("  A global proof may need factorization, Pell equations, or other methods.\n");
    end if;
end proc:


# ==========================================================
# STEP-BY-STEP MATHEMATICAL PROOF MODE (English, Student Level 2)
# Prints a proof of completeness ONLY when the exact divisor
# algorithm has established it. Otherwise gives an honest
# bounded-search explanation instead of a false proof.
# ==========================================================
ExplainProofMode := proc(a::integer,b::integer,c::integer,
                         d::integer,e::integer,f::integer,
                         info, proof, results::set,
                         complete::truefalse,
                         search_x::range, search_y::range,
                         max_solutions::posint)
    local P,D,K,A,B,C,E,F,G,k,h,N,Xsym,Ysym,Tsym,Usym,Vsym;
    local sol,xx,yy,Xv,Yv,Tv,Uv,Vv,idx,expected,divcount;
    local i,j,facpairs,u,v;
    P:=expand(a*x^2+b*x*y+c*y^2+d*x+e*y+f);
    D:=info["discriminant"];
    K:=info["determinant"];
    printf("\n========================================================\n");
    printf("  STEP-BY-STEP MATHEMATICAL PROOF MODE  (LEVEL 2)\n");
    printf("========================================================\n");
    printf("Question: Which integer pairs (x,y) satisfy P(x,y)=0?\n");
    print(P=0);
    printf("\nSTEP 1 / CLASSIFY THE EQUATION\n");
    printf("  What? Compute D=b^2-4*a*c and the conic determinant K.\n");
    printf("  Why? They describe the geometric type and degeneracy.\n");
    printf("  D=%a; K=%a.\n",D,K);
    printf("  Conclusion: %s\n",info["type"]);

    if not complete then
        printf("\nSTEP 2 / EXPLAIN THE AVAILABLE METHOD\n");
        printf("  An exact global divisor proof is not available for\n");
        printf("  this equation under the current algorithm/options.\n");
        if D>0 and not issqr(D) then
            printf("  D is not a square: a Pell-type method may be needed.\n");
        elif K=0 then
            printf("  The conic is degenerate: a separate line-factor proof may help.\n");
        else
            printf("  Other exact methods may be possible, but are not certified here.\n");
        end if;
        printf("\nSTEP 3 / CHECK THE RETURNED SOLUTIONS\n");
        idx:=0;
        for sol in results do
            xx:=eval(x,sol); yy:=eval(y,sol);
            idx:=idx+1;
            printf("  (%a,%a): P(x,y)=%a.\n",xx,yy,eval(P,sol));
        end do;
        printf("\nSTEP 4 / EXPLAIN WHAT IS (NOT) PROVED\n");
        printf("  Search rectangle: x=%a; y=%a.\n",search_x,search_y);
        printf("  At most %d solutions are returned.\n",max_solutions);
        printf("  Every returned point is verified, but the global\n");
        printf("  solution set is NOT proved complete.\n");
        if nops(results)>=max_solutions then
            printf("  The solution limit was reached; some points in the\n");
            printf("  search rectangle may also be omitted.\n");
        end if;
        return NULL;
    end if;

    if proof["swapped"] then
        A:=c; B:=b; C:=a; E:=e; F:=d; G:=f;
        printf("  To make the algebra work, set X=y and Y=x.\n");
    else
        A:=a; B:=b; C:=c; E:=d; F:=e; G:=f;
        printf("  Use X=x and Y=y.\n");
    end if;
    k:=proof["k"]; h:=proof["h"]; N:=proof["N"];
    printf("\nSTEP 2 / INTRODUCE AN EXACT INTEGER TRANSFORMATION\n");
    printf("  What? Set T=2*A*X+B*Y+E.\n");
    printf("  Why? Completing the square groups the quadratic terms.\n");
    printf("  Here A=%a, B=%a, C=%a, E=%a, F=%a, G=%a.\n",
              A,B,C,E,F,G);
    Tsym:=2*A*'X'+B*'Y'+E;
    Usym:=expand(k*Tsym-(D*'Y'+h));
    Vsym:=expand(k*Tsym+(D*'Y'+h));
    printf("  T = %a\n",Tsym);
    printf("  Let U=%a and V=%a.\n",Usym,Vsym);
    printf("  Here k=sqrt(D)=%a and h=B*E-2*A*F=%a.\n",k,h);
    printf("\nSTEP 3 / TURN THE EQUATION INTO A PRODUCT\n");
    printf("  Difference of squares: (p-q)*(p+q)=p^2-q^2.\n");
    printf("  For this equation the transformation yields U*V=%a.\n",N);
    printf("  Since X,Y are integers, U and V are integers.\n");
    printf("  The reverse formulas below also guarantee equivalence\n");
    printf("  after we retain only integer X,Y.\n");
    if a=1 and b=8 and c=7 and d=2 and e=11 and f=12 then
        printf("  A simpler, equivalent factorization in x,y is:\n");
        printf("    (2*x+2*y+3)*(2*x+14*y+1)=-45.\n");
        printf("  This is useful for understanding the divisor method.\n");
    end if;

    printf("\nSTEP 4 / WHY ARE THERE ONLY FINITELY MANY CASES?\n");
    divcount:=2*nops(numtheory:-divisors(abs(N)));
    printf("  Because U*V=%a and this is nonzero, U must divide %a.\n",N,N);
    printf("  There are %d signed possibilities for U.\n",divcount);
    printf("  For each U, the corresponding V is fixed: V=%a/U.\n",N);
    printf("  The algorithm examines EVERY such ordered factor pair.\n");

    printf("\nSTEP 5 / RECOVER THE ORIGINAL VARIABLES\n");
    printf("  T=(U+V)/(2*k) = (U+V)/%a.\n",2*k);
    printf("  Y=((V-U)/2-h)/D = ((V-U)/2-(%a))/%a.\n",h,D);
    printf("  X=(T-B*Y-E)/(2*A) = (T-(%a)*Y-(%a))/%a.\n",B,E,2*A);
    if proof["swapped"] then
        printf("  Convert back with x=Y, y=X.\n");
    else
        printf("  Here x=X, y=Y.\n");
    end if;
    printf("  Reject any factor pair producing non-integer x or y.\n");

    printf("\nSTEP 6 / SHOW THE SUCCESSFUL FACTOR PAIRS\n");
    idx:=0;
    for sol in results do
        xx:=eval(x,sol); yy:=eval(y,sol);
        if proof["swapped"] then Xv:=yy; Yv:=xx;
        else Xv:=xx; Yv:=yy; end if;
        Tv:=2*A*Xv+B*Yv+E;
        Uv:=k*Tv-(D*Yv+h);
        Vv:=k*Tv+(D*Yv+h);
        idx:=idx+1;
        printf("  %d) (x,y)=(%a,%a); (U,V)=(%a,%a); U*V=%a; P=%a.\n",
               idx,xx,yy,Uv,Vv,Uv*Vv,eval(P,sol));
    end do;
    if idx=0 then
        printf("  No factor pair produced integer x and y.\n");
    end if;
    printf("\nSTEP 7 / COMPLETENESS CONCLUSION\n");
    printf("  Every solution must correspond to a signed divisor U of %a.\n",N);
    printf("  Every possible signed divisor has been examined.\n");
    printf("  The inverse formulas were tested for integer coordinates,\n");
    printf("  and every accepted pair was checked in the original equation.\n");
    printf("  THEREFORE: exactly %d integer solutions exist globally.\n",idx);
    printf("  This is a proof, not merely a search within a graph window.\n");
    return NULL;
end proc:

# MAIN PROCEDURE
# The return value is the exact set when completeness is proved;
# otherwise it is a set bounded by search_x, search_y and max_solutions.
GeneralQuadraticSolver := proc(a::integer,b::integer,c::integer,
                               d::integer,e::integer,f::integer,
    {max_solutions::posint:=20,
     search_x::range:=-50..50, search_y::range:=-50..50,
     x_range::range:=-35..15, y_range::range:=-15..15,
     show_plot::truefalse:=true,
     show_details::truefalse:=true,
     show_general::truefalse:=true,
     show_education::truefalse:=true,
     education_level::string:="student",
     prove_complete::truefalse:=true,
     proof_mode::truefalse:=false})
    local P,info,D,K,proof,complete,results,sol,xx,yy,idx;
    local center,plot_result,general_list,verified;

    P:=expand(a*x^2+b*x*y+c*y^2+d*x+e*y+f);
    info:=ClassifyQuadratic(a,b,c,d,e,f);
    D:=info["discriminant"]; K:=info["determinant"];
    complete:=false;
    proof:=NULL;
    if prove_complete and D>0 and K<>0 and issqr(D) then
        proof:=SquareHyperbolaProof(a,b,c,d,e,f);
        if proof["proved"] then complete:=true; end if;
    end if;

    if complete then
        results:=proof["solutions"];
    else
        results:=FindIntegerSolutionsInBox(a,b,c,d,e,f,
                       search_x,search_y,max_solutions);
    end if;

    if show_details then
        printf("\n=======================================================\n");
        printf("       QUADRATIC DIOPHANTINE EQUATION SOLVER\n");
        printf("=======================================================\n");
        print(P=0);
        printf("Conic type         : %s\n",info["type"]);
        printf("Discriminant       : %a\n",D);
        printf("Conic determinant  : %a\n",K);
        printf("Integer solutions  : %d\n",nops(results));
        if complete then
            printf("Completeness       : PROVED - all integer solutions.\n");
        elif nops(results)>=max_solutions then
            printf("Completeness       : NOT PROVED - solution limit reached.\n");
        else
            printf("Completeness       : complete WITHIN the search box only.\n");
        end if;

        if show_education then
            if education_level = "student" then
                ExplainStudentLevel(a,b,c,d,e,f,info,proof,results,
                                    complete,search_x,search_y);
            elif education_level = "none" then
                # Skip the educational explanation.
            else
                error "education_level must be \"student\" or \"none\"";
            end if;
        end if;
        if proof_mode then
            ExplainProofMode(a,b,c,d,e,f,info,proof,results,
                             complete,search_x,search_y,max_solutions);
        end if;
        printf("\n--- VERIFIED INTEGER SOLUTIONS ---\n");
        idx:=0;
        for sol in results do
            xx:=eval(x,sol); yy:=eval(y,sol);
            verified:=eval(P,sol);
            idx:=idx+1;
            printf("  %d) x=%a, y=%a  [P(x,y)=%a]\n",
                   idx,xx,yy,verified);
        end do;
        if nops(results)=0 then printf("No solutions found by the selected method.\n"); end if;
        if not complete then
            printf("Search area: x=%a, y=%a; maximum returned=%d.\n",
                   search_x,search_y,max_solutions);
        end if;
    end if;

    if show_general then
        try
            general_list:=[isolve(P)];
            if show_details then
                printf("\n--- MAPLE ISOLVE OUTPUT ---\n");
                if nops(general_list)=0 then
                    printf("No parametric expression returned by isolve.\n");
                else
                    print(op(general_list));
                    printf("Any _Z-style parameter represents an integer.\n");
                end if;
            end if;
        catch:
            if show_details then
                printf("Maple isolve could not provide a result.\n");
            end if;
        end try;
    end if;

    if show_plot then
        plot_result:=PlotQuadraticFriendly(P,results,x_range,y_range,
                                          info["type"]);
        print(plot_result);
    end if;
    return results;
end proc:

# EXAMPLES (uncomment one at a time in Maple):
# GeneralQuadraticSolver(1,8,7,2,11,12);
# GeneralQuadraticSolver(1,8,7,2,11,12,x_range=-35..15);
# GeneralQuadraticSolver(1,0,-2,0,0,-1,show_plot=false);
# GeneralQuadraticSolver(1,0,-1,0,0,0);
# GeneralQuadraticSolver(1,0,0,0,-1,0);
# GeneralQuadraticSolver(1,0,1,0,0,-25);
# GeneralQuadraticSolver(1,0,1,0,0,1);
# GeneralQuadraticSolver(1,8,7,2,11,12,show_education=false);

# STUDENT LEVEL 2 EXAMPLE:
# GeneralQuadraticSolver(1,8,7,2,11,12,education_level="student",
#                        show_education=true,show_plot=true);

# STEP-BY-STEP PROOF MODE EXAMPLES:
# GeneralQuadraticSolver(1,8,7,2,11,12,
#      proof_mode=true,show_education=false,show_plot=true);
# GeneralQuadraticSolver(1,0,-2,0,0,-1,
#      proof_mode=true,show_education=false,show_plot=false);
# GeneralQuadraticSolver(1,0,1,0,0,-25,
#      proof_mode=true,show_education=false,show_plot=false);


# ==========================================================
# INTERACTIVE MAPLET FRONT END
# Level 2 English; no solver changes required.
# ==========================================================

GQS_M_ReadInt := proc(id)
    local txt;
    txt := StringTools:-Trim(convert(Maplets:-Tools:-Get(id),string));
    if not StringTools:-RegMatch("^-?[0-9]+$",txt) then
        error "Input %1 must be an integer (received %2)",id,txt;
    end if;
    return parse(txt);
end proc:

GQS_M_ReadBool := proc(id)
    local z;
    z:=Maplets:-Tools:-Get(id);
    return evalb(z=true or z="true");
end proc:

GQS_M_Inputs := proc()
    local a,b,c,d,e,f,xlo,xhi,ylo,yhi,pxlo,pxhi,pylo,pyhi,limit;
    a:=GQS_M_ReadInt('coefA'); b:=GQS_M_ReadInt('coefB');
    c:=GQS_M_ReadInt('coefC'); d:=GQS_M_ReadInt('coefD');
    e:=GQS_M_ReadInt('coefE'); f:=GQS_M_ReadInt('coefF');
    xlo:=GQS_M_ReadInt('searchXmin'); xhi:=GQS_M_ReadInt('searchXmax');
    ylo:=GQS_M_ReadInt('searchYmin'); yhi:=GQS_M_ReadInt('searchYmax');
    pxlo:=GQS_M_ReadInt('plotXmin'); pxhi:=GQS_M_ReadInt('plotXmax');
    pylo:=GQS_M_ReadInt('plotYmin'); pyhi:=GQS_M_ReadInt('plotYmax');
    limit:=GQS_M_ReadInt('solLimit');
    if xlo>xhi or ylo>yhi or pxlo>=pxhi or pylo>=pyhi then
        error "Check ranges: minimum must not exceed maximum";
    end if;
    if limit<=0 then error "Solution limit must be positive"; end if;
    return [a,b,c,d,e,f,xlo..xhi,ylo..yhi,pxlo..pxhi,pylo..pyhi,limit];
end proc:

GQS_M_Report := proc()
    local z,a,b,c,d,e,f,sx,sy,px,py,lim,info,proof,complete,sols,P,txt,i,sol;
    local D,k,h,N,A,B,C,E,F,G,invX,invY,xx,yy,den0,cx0,cy0,val0,kind0;
    try
        z:=GQS_M_Inputs();
        a:=z[1]; b:=z[2]; c:=z[3]; d:=z[4]; e:=z[5]; f:=z[6];
        sx:=z[7]; sy:=z[8]; px:=z[9]; py:=z[10]; lim:=z[11];
        P:=expand(a*x^2+b*x*y+c*y^2+d*x+e*y+f);
        info:=ClassifyQuadratic(a,b,c,d,e,f);
        D:=info["discriminant"];
        proof:=NULL; complete:=false;
        if GQS_M_ReadBool('proofCheck') and D>0
           and info["determinant"]<>0 and issqr(D) then
            proof:=SquareHyperbolaProof(a,b,c,d,e,f);
            complete:=evalb(proof["proved"]=true);
        end if;
        if complete then
            sols:=proof["solutions"];
        else
            sols:=FindIntegerSolutionsInBox(a,b,c,d,e,f,sx,sy,lim);
        end if;
        txt:=sprintf("DOMAIN: INTEGERS (Z) - DIOPHANTINE SOLUTIONS\\n");
        txt:=cat(txt,"The real conic and complex curve use the same equation, but different number systems.\\n\\n");
        txt:=cat(txt,sprintf("Equation: %a = 0\n",P));
        txt:=cat(txt,sprintf("Classification: %s\n",info["type"]));
        txt:=cat(txt,sprintf("Discriminant: %a   Determinant: %a\n\n",
                             D,info["determinant"]));
        txt:=cat(txt,"------------------------------\nSTEP 1 - IDENTIFY THE CURVE\n------------------------------\n");
        txt:=cat(txt,"The discriminant tells us which family of curves we have.\n");
        txt:=cat(txt,"The determinant tells us if the curve splits into simpler parts.\n\n");

        if info["type"]="IMAGINARY ELLIPSE (NO REAL POINTS)" then
            den0:=4*a*c-b^2;
            cx0:=(b*e-2*c*d)/den0;
            cy0:=(b*d-2*a*e)/den0;
            val0:=simplify(eval(P,{x=cx0,y=cy0}));
            if a>0 then kind0:="minimum"; else kind0:="maximum"; end if;
            txt:=cat(txt,"------------------------------\n");
            txt:=cat(txt,"STEP 2 - FIND THE CLOSEST VALUE TO ZERO\n");
            txt:=cat(txt,"------------------------------\n");
            txt:=cat(txt,sprintf("The expression has its %s at:\n",kind0));
            txt:=cat(txt,sprintf("  x = %a, y = %a\n",cx0,cy0));
            txt:=cat(txt,sprintf("  P = %a (approximately %.4f)\n\n",
                                 val0,evalf(val0)));
            txt:=cat(txt,"------------------------------\n");
            txt:=cat(txt,"STEP 3 - WHY IS THERE NO SOLUTION?\n");
            txt:=cat(txt,"------------------------------\n");
            if a>0 then
                txt:=cat(txt,"This is the smallest possible value of P.\n");
                txt:=cat(txt,"It is greater than zero, so P can never equal zero.\n\n");
            else
                txt:=cat(txt,"This is the largest possible value of P.\n");
                txt:=cat(txt,"It is less than zero, so P can never equal zero.\n\n");
            end if;
            txt:=cat(txt,"------------------------------\n");
            txt:=cat(txt,"STEP 4 - UNDERSTAND THE GRAPH\n");
            txt:=cat(txt,"------------------------------\n");
            txt:=cat(txt,"Think of contour lines on a landscape map.\n");
            txt:=cat(txt,"Each ellipse joins points with the SAME value of P.\n");
            txt:=cat(txt,"The blue curve is closest to the red extremum.\n");
            txt:=cat(txt,"The green and purple curves show larger distances.\n");
            txt:=cat(txt,"IMPORTANT: these curves are NOT solutions of P=0.\n");
            txt:=cat(txt,"The red dot marks the lowest or highest value of P.\n\n");
            txt:=cat(txt,"CONCLUSION: NO REAL OR INTEGER SOLUTIONS EXIST.\n");
            return txt;
        end if;

        if complete then
            txt:=cat(txt,"------------------------------\nSTEP 2 - REWRITE THE EQUATION\n------------------------------\n");
            if proof["swapped"] then
                A:=c; B:=b; C:=a; E:=e; F:=d; G:=f;
                txt:=cat(txt,"For this transformation, X=y and Y=x.\n");
            else
                A:=a; B:=b; C:=c; E:=d; F:=e; G:=f;
                txt:=cat(txt,"For this transformation, X=x and Y=y.\n");
            end if;
            k:=proof["k"]; h:=proof["h"]; N:=proof["N"];
            txt:=cat(txt,sprintf("Define T = %a*X + %a*Y + %a.\n",2*A,B,E));
            txt:=cat(txt,sprintf("D=%a, k=sqrt(D)=%a, h=%a.\n",D,k,h));
            txt:=cat(txt,sprintf("Then (k*T-D*Y-h)(k*T+D*Y+h) = %a.\n",N));
            txt:=cat(txt,"------------------------------\nSTEP 3 - CHECK THE POSSIBLE FACTORS\n------------------------------\n");
            txt:=cat(txt,sprintf("Both factors are integers. They must multiply to %a.\n",N));
            txt:=cat(txt,"Only finitely many signed divisor pairs are possible.\n");
            txt:=cat(txt,"------------------------------\nSTEP 4 - FIND x AND y\n------------------------------\n");
            txt:=cat(txt,"For each pair (u,v), compute:\n");
            txt:=cat(txt,sprintf("  T=(u+v)/(%a)\n",2*k));
            txt:=cat(txt,sprintf("  Y=((v-u)/2-(%a))/%a\n",h,D));
            txt:=cat(txt,sprintf("  X=(T-(%a)*Y-(%a))/%a\n",B,E,2*A));
            txt:=cat(txt,"Keep only pairs that give integer coordinates.\n\n");
        else
            txt:=cat(txt,"------------------------------\nSTEP 2 - SEARCH FOR INTEGER PAIRS\n------------------------------\n");
            txt:=cat(txt,sprintf("Checking x=%a and y=%a, at most %d solutions.\n",sx,sy,lim));
            txt:=cat(txt,"An unrestricted completeness proof is not available with the chosen method.\n\n");
        end if;
        txt:=cat(txt,"------------------------------\nSTEP 5 - CHECK EACH SOLUTION\n------------------------------\n");
        i:=0;
        for sol in sols do
            i:=i+1;
            xx:=eval(x,sol); yy:=eval(y,sol);
            txt:=cat(txt,sprintf("%d) x=%a, y=%a; P(x,y)=%a\n",i,xx,yy,eval(P,sol)));
        end do;
        if i=0 then txt:=cat(txt,"No integer solutions returned.\n"); end if;
        txt:=cat(txt,"\n------------------------------\nSTEP 6 - WHAT HAVE WE PROVED?\n------------------------------\n");
        if complete then
            txt:=cat(txt,sprintf("GLOBAL PROOF: Exactly %d integer solutions exist.\n",i));
            txt:=cat(txt,"Every possible factor pair was checked, not just points in the search box.\n");
        elif i>=lim then
            txt:=cat(txt,"ONLY A LIMITED SEARCH: solution cap reached; other pairs may exist even inside the box.\n");
        else
            txt:=cat(txt,"SEARCH-BOX COMPLETE: all integer pairs in this box were checked.\n");
            txt:=cat(txt,"Additional solutions outside the box may exist.\n");
        end if;
        return txt;
    catch:
        return cat("INPUT OR COMPUTATION ERROR:\n",lasterror,
                   "\n\nCheck the integer entries and search ranges.");
    end try;
end proc:

# Send true separate lines to the multi-line TextBox.
GQS_M_ShowReport := proc()
    local txt, lines, line;
    txt := GQS_M_Report();
    Maplets:-Tools:-Set('ResultBox'(value) = "");
    lines := StringTools:-Split(txt, sprintf("%c", 10));
    for line in lines do
        Maplets:-Tools:-Set('ResultBox'(appendline) = line);
    end do;
    return NULL;
end proc:

GQS_M_Plot := proc()
    local z,a,b,c,d,e,f,info,D,proof,sols,P;
    z:=GQS_M_Inputs();
    a:=z[1]; b:=z[2]; c:=z[3]; d:=z[4]; e:=z[5]; f:=z[6];
    P:=expand(a*x^2+b*x*y+c*y^2+d*x+e*y+f);
    info:=ClassifyQuadratic(a,b,c,d,e,f);
    D:=info["discriminant"];
    if GQS_M_ReadBool('proofCheck') and D>0
       and info["determinant"]<>0 and issqr(D) then
        proof:=SquareHyperbolaProof(a,b,c,d,e,f);
        if proof["proved"] then
            sols:=proof["solutions"];
        else
            sols:=FindIntegerSolutionsInBox(a,b,c,d,e,f,z[7],z[8],z[11]);
        end if;
    else
        sols:=FindIntegerSolutionsInBox(a,b,c,d,e,f,z[7],z[8],z[11]);
    end if;
    return PlotQuadraticFriendly(P,sols,z[9],z[10],info["type"]);
end proc:

# BEGINNER-FRIENDLY INPUT: example presets and a short help guide.
GQS_M_LoadCoefficients := proc(values)
    local ids, i;
    ids := ['coefA','coefB','coefC','coefD','coefE','coefF'];
    for i from 1 to 6 do
        Maplets:-Tools:-Set(ids[i](value) = convert(values[i],string));
    end do;
    return NULL;
end proc:

GQS_M_ExampleHyperbola := proc()
    GQS_M_LoadCoefficients([1,8,7,2,11,12]);
end proc:
GQS_M_ExampleEllipse := proc()
    GQS_M_LoadCoefficients([1,0,1,0,0,-25]);
end proc:
GQS_M_ExampleParabola := proc()
    GQS_M_LoadCoefficients([1,0,0,0,-1,0]);
end proc:
GQS_M_ExamplePell := proc()
    GQS_M_LoadCoefficients([1,0,-2,0,0,-1]);
end proc:
GQS_M_ExampleTwoLines := proc()
    GQS_M_LoadCoefficients([1,0,-1,0,0,0]);
end proc:

GQS_M_InputHelp := proc()
    local lines, line;
    lines := [
      "HOW TO ENTER AN EQUATION",
      "",
      "Easiest option: click one of the example buttons at the top.",
      "The six boxes then fill in automatically.",
      "",
      "GENERAL FORM: a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0",
      "",
      "a = coefficient before x^2",
      "b = coefficient before x*y",
      "c = coefficient before y^2",
      "d = coefficient before x",
      "e = coefficient before y",
      "f = constant term (no x or y)",
      "",
      "Missing terms? Enter 0 for those coefficients.",
      "",
      "EXAMPLE: x^2 + y^2 = 25",
      "Move everything to the left: x^2 + y^2 - 25 = 0.",
      "Enter: a=1, b=0, c=1, d=0, e=0, f=-25.",
      "",
      "Next: click SOLVE + EXPLAIN, then DRAW GRAPH.",
      "The program searches only for INTEGER solutions."
    ];
    Maplets:-Tools:-Set('ResultBox'(value) = "");
    for line in lines do
       Maplets:-Tools:-Set('ResultBox'(appendline) = line);
    end do;
    return NULL;
end proc:


# ------------------------------------------------------------
# OPTIONAL COMPLEX-NUMBER EXPLORATION (a two-dimensional SLICE).
# Fix y=0 when possible. The displayed complex roots belong to
# the selected slice; they do NOT constitute the whole complex conic.
# ------------------------------------------------------------
GQS_M_ComplexData := proc()
    local z,a,b,c,d,e,f,P,eq,var,fixed,poly,roots;
    z := GQS_M_Inputs();
    a:=z[1]; b:=z[2]; c:=z[3]; d:=z[4]; e:=z[5]; f:=z[6];
    P:=expand(a*x^2+b*x*y+c*y^2+d*x+e*y+f);
    if a<>0 or d<>0 then
        fixed := "y = 0";
        var := x;
        poly := expand(subs(y=0,P));
    elif c<>0 or e<>0 then
        fixed := "x = 0";
        var := y;
        poly := expand(subs(x=0,P));
    else
        # On the coordinate axes P reduces to f.
        # If f=0, both axes are solutions, not isolated roots.
        return [P, "both axes", x, f, []];
    end if;
    if degree(poly,var)=0 then
        roots:=[];
    else
        roots := [solve(poly=0,var)];
    end if;
    return [P,fixed,var,poly,roots];
end proc:

GQS_M_ComplexReport := proc()
    local data,P,fixed,var,poly,roots,lines,line,r,v;
    data:=GQS_M_ComplexData();
    P:=data[1]; fixed:=data[2]; var:=data[3];
    poly:=data[4]; roots:=data[5];
    lines := [
        "DOMAIN: COMPLEX NUMBERS (C) - ONE SLICE",
        "These roots are NOT a plot of the full complex conic.",
        "A real root is also complex (its imaginary part is 0).",
        "",
        sprintf("Original equation: %a = 0",P),
        "",
        sprintf("STEP 1 - Choose the slice %s.",fixed),
        "We keep one variable equal to zero.",
        "This leaves a polynomial in just one variable.",
        sprintf("Slice equation: %a = 0",poly),
        "",
        "STEP 2 - Solve over the complex numbers"
    ];
    if degree(poly,var)=0 then
        if poly=0 then
            lines:=[op(lines),
              "Every point on the chosen coordinate axis is a solution.",
              "There is no finite set of isolated roots to plot."];
        else
            lines:=[op(lines),
              "This coordinate-axis slice has no solutions.",
              "Other complex points may still exist away from this slice."];
        end if;
    elif nops(roots)=0 then
        lines:=[op(lines),"No roots returned for this slice."];
    else
        for r in roots do
            lines:=[op(lines),sprintf("%a = %a",var,r),
                    sprintf("Approximation: %.5f %+.5f*i",
                       evalf(Re(r)),evalf(Im(r)))];
        end do;
    end if;
    lines:=[op(lines),"",
        "STEP 3 - INTERPRET THE GRAPH",
        "Horizontal coordinate: real part.",
        "Vertical coordinate: imaginary part.",
        "The points are roots in the complex plane.",
        "IMPORTANT: this is a slice, NOT the original x-y plot.",
        "It is not a drawing of the entire complex conic."
    ];
    Maplets:-Tools:-Set('ResultBox'(value)="");
    for line in lines do
        Maplets:-Tools:-Set('ResultBox'(appendline)=line);
    end do;
    return NULL;
end proc:

GQS_M_ComplexPlot := proc()
    local data, roots, r, pts, rr, ii, m, ax;
    data:=GQS_M_ComplexData();
    roots:=data[5];
    pts:=[];
    for r in roots do
        rr:=evalf(Re(r)); ii:=evalf(Im(r));
        if type(rr,numeric) and type(ii,numeric) then
            pts:=[op(pts),[rr,ii]];
        end if;
    end do;
    if nops(pts)=0 then
        return plots:-display(plot([], -2..2),
          title="No isolated roots in this slice",
          labels=["Real part","Imaginary part"]);
    end if;
    m:=max(1,op(map(u->max(abs(u[1]),abs(u[2])),pts)));
    ax:=evalf(1.3*m);
    return plots:-display(
       plots:-pointplot(pts,color=red,symbol=solidcircle,
                       symbolsize=17,legend="Complex roots"),
       axes=normal,view=[-ax..ax,-ax..ax],
       scaling=constrained,
       labels=["Real part","Imaginary part"],
       title="COMPLEX ROOTS - coordinate-axis slice");
end proc:


# ------------------------------------------------------------
# COMPLEX ALGEBRAIC CURVE: 3D REAL PROJECTION
#
# Independent complex coordinate z=u+I*v.
# The dependent complex coordinate is solved quadratically.
# 3D plot: (Re(z), Im(z), Re(w)) for BOTH square-root branches.
# NOTE: Im(w) is omitted: this is a projection from R^4 to R^3.
# ------------------------------------------------------------
GQS_M_ComplexCurveData := proc()
    local z,a,b,c,d,e,f,ind,dep,A0,B0,C0,w;
    z:=GQS_M_Inputs();
    a:=z[1]; b:=z[2]; c:=z[3]; d:=z[4]; e:=z[5]; f:=z[6];
    if c<>0 then
        ind:=x; dep:=y;
        A0:=c;
        B0:=b*Z+e;
        C0:=a*Z^2+d*Z+f;
    elif a<>0 then
        ind:=y; dep:=x;
        A0:=a;
        B0:=b*Z+d;
        C0:=c*Z^2+e*Z+f;
    else
        return NULL;
    end if;
    return [ind,dep,A0,B0,C0];
end proc:

GQS_M_ComplexCurveReport := proc()
    local data,lines,line,z,a,b,c,d,e,f,P;
    z:=GQS_M_Inputs();
    a:=z[1]; b:=z[2]; c:=z[3];
    d:=z[4]; e:=z[5]; f:=z[6];
    P:=expand(a*x^2+b*x*y+c*y^2+d*x+e*y+f);
    data:=GQS_M_ComplexCurveData();
    if [a,b,c,d,e,f]=[1,0,-1,0,0,0] then
        lines:=[
          "COMPLEX ALGEBRAIC CURVE - TWO LINES",
          "These stay two complex lines, NOT a hyperbola.",
          "In this projection they appear as two planes.",
          "Equation: x^2 - y^2 = 0",
          "",
          "STEP 1 - FACTORIZE",
          "(x-y)*(x+y)=0",
          "",
          "STEP 2 - REAL GEOMETRY",
          "Line 1: y=x; Line 2: y=-x.",
          "There are INFINITELY MANY integer solutions:",
          "(x,y)=(n,n) or (n,-n), for every integer n.",
          "",
          "STEP 3 - EXTEND TO COMPLEX NUMBERS",
          "Write z=x and w=y with complex coordinates.",
          "Then w=z or w=-z defines two complex lines.",
          "",
          "STEP 4 - UNDERSTAND THE 3D PROJECTION",
          "CYAN: projected complex plane w=z.",
          "MAGENTA: projected complex plane w=-z.",
          "RED: the two original real lines (Im(z)=Im(w)=0).",
          "BLACK: some integer points in the selected bounds.",
          "The complete complex curve lives in four real dimensions.",
          "This projection omits Im(w).",
          "These are still two complex lines (a degenerate conic),",
          "although the 3D projections appear as planes."
        ];
    elif data=NULL then
        lines:=[
          "DOMAIN: COMPLEX NUMBERS (C)",
          "This is a quadratic algebraic curve, not an",
          "ordinary Euclidean ellipse or hyperbola.",
          "Pell refers only to INTEGER solutions.",
          "",
          "This example has no squared variable.",
          "The quadratic two-sheet display needs a nonzero",
          "coefficient of x^2 or y^2.",
          "Use EXPLORE COMPLEX SOLUTIONS instead."
        ];
    else
        lines:=[
          "DOMAIN: COMPLEX NUMBERS (C)",
          "This is a quadratic algebraic curve, not an",
          "ordinary Euclidean ellipse or hyperbola.",
          "Pell refers only to INTEGER solutions.",
          "",
          "STEP 1 - Allow both variables to be complex.",
          sprintf("Let %a = u + I*v.",data[1]),
          "",
          sprintf("STEP 2 - Solve the quadratic equation for %a.",data[2]),
          "There are two square-root branches in general.",
          "",
          "STEP 3 - Display both branches in three dimensions.",
          sprintf("Axis 1: Re(%a) = u",data[1]),
          sprintf("Axis 2: Im(%a) = v",data[1]),
          sprintf("Axis 3: Re(%a)",data[2]),
          "",
          "The omitted coordinate is the imaginary part",
          "of the dependent variable.",
          "",
          "IMPORTANT: this is a 3D PROJECTION of a complex",
          "curve in four real dimensions, not the full object.",
          "The two branches can overlap in this projection.",
          "",
          "IMPORTANT MATHEMATICAL DISTINCTION",
          "The real graph is a conic in the ordinary x-y plane.",
          "Over C, the same polynomial defines a complex",
          "algebraic curve of degree two.",
          "The displayed 3D projection is NOT a Euclidean conic.",
          "A Pell equation concerns integer points, not",
          "the entire complex solution set.",
          "",
          "RED curve: original real conic (where both imaginary parts=0).",
          "BLACK dots: integer solutions inside the displayed domain.",
          "CYAN/MAGENTA: projected complex branches.",
          "Try dragging the 3D graph to rotate it."
        ];
    end if;
    Maplets:-Tools:-Set('ResultBox'(value)="");
    for line in lines do
        Maplets:-Tools:-Set('ResultBox'(appendline)=line);
    end do;
    return NULL;
end proc:

GQS_M_ComplexCurvePlot := proc()
    local data,A0,B0,C0,wplus,wminus,disc,zexpr,u,v,qplus,qminus,
          p1,p2,layers,zz,dd,ww1,ww2,branch1,branch2,j,tmin,tmax,
          dval,w1val,w2val,z,aa,bb,cc,ee,ff,gg,P,sols,sol,xx,yy,
          segments1,segments2,
          zval,wval,intpts,xyRange,iyRange,limit,ind,plotCaption;
    data:=GQS_M_ComplexCurveData();
    if data=NULL then
        return plots:-display(plot([], -6..6),
               title="3D projection unavailable for this equation");
    end if;
    A0:=data[3]; B0:=data[4]; C0:=data[5];
    ind:=data[1];
    tmin:=-10; tmax:=10;
    zexpr:=u+I*v;
    disc:=expand(B0^2-4*A0*C0);
    wplus := subs(Z=zexpr,(-B0+sqrt(disc))/(2*A0));
    wminus:= subs(Z=zexpr,(-B0-sqrt(disc))/(2*A0));

    qplus := evalc(Re(wplus));
    qminus:= evalc(Re(wminus));
    p1:=plot3d([u,v,qplus],u=tmin..tmax,v=tmin..tmax,
              grid=[35,35],color=cyan,
              style=surface,transparency=0.35);
    p2:=plot3d([u,v,qminus],u=tmin..tmax,v=tmin..tmax,
              grid=[35,35],color=magenta,
              style=surface,transparency=0.55);
    layers:=[p1,p2];

    # Real points have Im(z)=Im(w)=0. For each real z with
    # nonnegative discriminant, plot both real quadratic roots.
    branch1:=[]; branch2:=[];
    segments1:=[]; segments2:=[];
    for j from 0 to 600 do
        zz:=tmin+(tmax-tmin)*j/600;
        dval:=eval(disc,Z=zz);
        if evalb(dval>=0) then
            w1val:=evalf(subs(Z=zz,(-B0+sqrt(disc))/(2*A0)));
            w2val:=evalf(subs(Z=zz,(-B0-sqrt(disc))/(2*A0)));
            if type(w1val,numeric) and type(w2val,numeric) then
                branch1:=[op(branch1),[evalf(zz),0,w1val]];
                branch2:=[op(branch2),[evalf(zz),0,w2val]];
            end if;
        else
            if nops(branch1)>1 then
                segments1:=[op(segments1),branch1];
                segments2:=[op(segments2),branch2];
            end if;
            branch1:=[]; branch2:=[];
        end if;
    end do;
    if nops(branch1)>1 then
        segments1:=[op(segments1),branch1];
        segments2:=[op(segments2),branch2];
    end if;
    for sol in segments1 do
        layers:=[op(layers),
                plots:-spacecurve(sol,color=red,thickness=4)];
    end do;
    for sol in segments2 do
        layers:=[op(layers),
                plots:-spacecurve(sol,color=red,thickness=4)];
    end do;

    # Black markers: integer points within BOTH the search box and
    # the real-axis domain of the displayed 3D projection.
    z:=GQS_M_Inputs();
    aa:=z[1]; bb:=z[2]; cc:=z[3];
    ee:=z[4]; ff:=z[5]; gg:=z[6];
    P:=aa*x^2+bb*x*y+cc*y^2+ee*x+ff*y+gg;
    sols:=FindIntegerSolutionsInBox(aa,bb,cc,ee,ff,gg,z[7],z[8],z[11]);
    intpts:=[];
    for sol in sols do
        xx:=eval(x,sol); yy:=eval(y,sol);
        if ind=x then zval:=xx; wval:=yy;
        else zval:=yy; wval:=xx; end if;
        if zval>=tmin and zval<=tmax then
            intpts:=[op(intpts),[zval,0,wval]];
        end if;
    end do;
    if nops(intpts)>0 then
        layers:=[op(layers),
           plots:-pointplot3d(intpts,color=black,
                        symbol=solidcircle,symbolsize=14)];
    end if;
    if [aa,bb,cc,ee,ff,gg]=[1,0,-1,0,0,0] then
        plotCaption := "CYAN: w=z  |  MAGENTA: w=-z  |  RED: y=x and y=-x  |  BLACK: sample integer points\nThese remain TWO COMPLEX LINES (shown as planes), NOT a hyperbola. Integer solutions: (n,n) or (n,-n). Im(w) omitted.";
    else
        plotCaption :=  "CYAN/MAGENTA: complex branches  |  RED: real conic  |  BLACK: integer points\nThis is NOT an ordinary conic or a Pell graph. The complex curve has 4 real dimensions; the 3D projection omits Im(w).";
    end if;
    return plots:-display(layers,axes=boxed,
             labels=["Re(z)","Im(z)","Re(w)"],
             title=sprintf("OVER C: COMPLEX QUADRATIC CURVE  |  %a = 0",P),
             caption=plotCaption,
             orientation=[-120,65]);
end proc:


# EDUCATIONAL BRIDGE: explains why the same equation is studied over Z, R and C.
# This is explanation, not a claim that the complex visualization proves anything
# about integer completeness.
GQS_M_MathBridge := proc()
    local z,a,b,c,d,e,f,P,info,D,kind,lines,line;
    z:=GQS_M_Inputs();
    a:=z[1]; b:=z[2]; c:=z[3];
    d:=z[4]; e:=z[5]; f:=z[6];
    P:=expand(a*x^2+b*x*y+c*y^2+d*x+e*y+f);
    info:=ClassifyQuadratic(a,b,c,d,e,f);
    D:=info["discriminant"];
    kind:=info["type"];
    lines:=[
      "ONE EQUATION - THREE MATHEMATICAL WORLDS",
      "",
      sprintf("Equation: %a = 0",P),
      "",
      "1) INTEGERS (Z) - NUMBER THEORY",
      "The Diophantine problem asks for x and y in Z.",
      "Its answers are individual lattice points.",
      "A proof of ALL integer solutions requires number theory.",
      "",
      "2) REALS (R) - ANALYTIC GEOMETRY",
      sprintf("Real classification: %s",kind),
      sprintf("Discriminant: %a",D),
      "The real graph shows all real points of P(x,y)=0.",
      "The integer points lie on this graph, if they exist.",
      "",
      "3) COMPLEX NUMBERS (C) - ALGEBRAIC GEOMETRY",
      "Now x and y may each have real and imaginary parts.",
      "The equation is unchanged; its solution domain expands.",
      "The complex curve has TWO real dimensions in FOUR",
      "real coordinate dimensions.",
      "Our 3D plot is a projection, NOT an ordinary conic.",
      "",
      "WHAT DOES THIS ADD?",
      "The complex plot does NOT find more integer solutions.",
      "It explains the larger algebraic object containing",
      "the real curve (when real points exist)."
    ];
    if [a,b,c,d,e,f]=[1,0,-2,0,0,-1] then
        lines:=[op(lines),"",
          "PELL CONNECTION",
          "x^2-2*y^2=1 also means",
          "(x+sqrt(2)*y)*(x-sqrt(2)*y)=1.",
          "Integer solutions are governed by units in Z[sqrt(2)].",
          "The complex plot visualizes the algebraic curve,",
          "but does not replace this integer argument."];
    elif [a,b,c,d,e,f]=[1,0,-1,0,0,0] then
        lines:=[op(lines),"",
          "TWO-LINES CONNECTION",
          "(x-y)*(x+y)=0 gives y=x or y=-x.",
          "Over C the same factorization gives two complex lines.",
          "These appear as planes in the 3D projection.",
          "Integer solutions are (n,n) and (n,-n), n in Z."];
    elif kind="IMAGINARY ELLIPSE (NO REAL POINTS)" then
        lines:=[op(lines),"",
          "IMAGINARY-ELLIPSE CONNECTION",
          "P(x,y)=0 has no real points, hence no integer points.",
          "It may still have points over C.",
          "Colored real contour lines are NOT complex solutions."];
    elif kind="REAL ELLIPSE (NON-DEGENERATE)" then
        lines:=[op(lines),"",
          "ELLIPSE CONNECTION",
          "The real ellipse is the subset with both imaginary",
          "parts zero. Complex extensions need not look elliptical",
          "after projection."];
    end if;
    lines:=[op(lines),"",
        "STUDY TIP",
        "First select INTEGER SOLUTIONS + PROOF.",
        "Next select REAL CONIC - 2D.",
        "Finally try the ADVANCED complex exploration."];
    Maplets:-Tools:-Set('ResultBox'(value)="");
    for line in lines do
        Maplets:-Tools:-Set('ResultBox'(appendline)=line);
    end do;
    return NULL;
end proc:

GQS_Maplet := proc()
    local M;
    M:=Maplets:-Elements:-Maplet(
        'onstartup'='startupAction',
        Maplets:-Elements:-Window['MainWindow'](
            'title'="Quadratic Diophantine Tutor - Step-by-Step Proof (English)",
            'width'=1260, 'height'=925,
            Maplets:-Elements:-BoxColumn(
                'inset'=8, 'spacing'=8,
                Maplets:-Elements:-Label("ONE EQUATION - THREE MATHEMATICAL WORLDS"),
                Maplets:-Elements:-Label("INTEGERS (Z): Diophantine points  |  REALS (R): ordinary conic  |  COMPLEX (C): algebraic curve, not a Euclidean conic"),
                Maplets:-Elements:-Label("A Pell equation describes an integer problem; the complex 3D surface is NOT a Pell plot."),
                Maplets:-Elements:-Label("STEP A: Choose an example (the six numbers fill in automatically)"),
                Maplets:-Elements:-BoxRow(
                    Maplets:-Elements:-Button("HYPERBOLA", 'onclick'='exampleHyperbolaAction'),
                    Maplets:-Elements:-Button("ELLIPSE", 'onclick'='exampleEllipseAction'),
                    Maplets:-Elements:-Button("PARABOLA", 'onclick'='exampleParabolaAction'),
                    Maplets:-Elements:-Button("PELL", 'onclick'='examplePellAction'),
                    Maplets:-Elements:-Button("TWO LINES", 'onclick'='exampleLinesAction'),
                    Maplets:-Elements:-Button("HOW TO ENTER AN EQUATION?", 'onclick'='inputHelpAction')
                ),
                Maplets:-Elements:-Label("STEP B (optional): Edit a,b,c,d,e,f for a*x^2+b*x*y+c*y^2+d*x+e*y+f=0"),
                Maplets:-Elements:-Label("a = x^2   |   b = x*y   |   c = y^2   |   d = x   |   e = y   |   f = constant   |   Missing term = 0"),
                Maplets:-Elements:-BoxRow(
                    Maplets:-Elements:-Label("a"),
                    Maplets:-Elements:-TextField['coefA'](value="1",width=5),
                    Maplets:-Elements:-Label("b"),
                    Maplets:-Elements:-TextField['coefB'](value="8",width=5),
                    Maplets:-Elements:-Label("c"),
                    Maplets:-Elements:-TextField['coefC'](value="7",width=5),
                    Maplets:-Elements:-Label("d"),
                    Maplets:-Elements:-TextField['coefD'](value="2",width=5),
                    Maplets:-Elements:-Label("e"),
                    Maplets:-Elements:-TextField['coefE'](value="11",width=5),
                    Maplets:-Elements:-Label("f"),
                    Maplets:-Elements:-TextField['coefF'](value="12",width=5)
                ),
                Maplets:-Elements:-BoxRow(
                    Maplets:-Elements:-Label("Search x:"),
                    Maplets:-Elements:-TextField['searchXmin'](value="-50",width=5),
                    Maplets:-Elements:-Label("to"),
                    Maplets:-Elements:-TextField['searchXmax'](value="50",width=5),
                    Maplets:-Elements:-Label("Search y:"),
                    Maplets:-Elements:-TextField['searchYmin'](value="-50",width=5),
                    Maplets:-Elements:-Label("to"),
                    Maplets:-Elements:-TextField['searchYmax'](value="50",width=5),
                    Maplets:-Elements:-Label("Max solutions:"),
                    Maplets:-Elements:-TextField['solLimit'](value="20",width=5)
                ),
                Maplets:-Elements:-BoxRow(
                    Maplets:-Elements:-Label("Plot x:"),
                    Maplets:-Elements:-TextField['plotXmin'](value="-35",width=5),
                    Maplets:-Elements:-Label("to"),
                    Maplets:-Elements:-TextField['plotXmax'](value="15",width=5),
                    Maplets:-Elements:-Label("Plot y:"),
                    Maplets:-Elements:-TextField['plotYmin'](value="-15",width=5),
                    Maplets:-Elements:-Label("to"),
                    Maplets:-Elements:-TextField['plotYmax'](value="15",width=5),
                    Maplets:-Elements:-CheckBox['proofCheck']("Try exact completeness proof",value=true)
                ),
                Maplets:-Elements:-Label("CORE LEARNING - START HERE"),
                Maplets:-Elements:-BoxRow(
                    Maplets:-Elements:-Button("1. INTEGER SOLUTIONS + PROOF", 'onclick'='runAction'),
                    Maplets:-Elements:-Button("2. REAL CONIC - 2D", 'onclick'='plotAction'),
                    Maplets:-Elements:-Button("3. EXPLAIN THE CONNECTION", 'onclick'='bridgeAction')
                ),
                Maplets:-Elements:-Label("ADVANCED EXPLORATION - OPTIONAL (complex numbers)"),
                Maplets:-Elements:-BoxRow(
                    Maplets:-Elements:-Button("COMPLEX ROOTS - SLICE",
                                             'onclick'='complexAction'),
                    Maplets:-Elements:-Button("COMPLEX CURVE - 3D",
                                             'onclick'='complexCurveAction'),
                    Maplets:-Elements:-Button("CLOSE", 'onclick'=Maplets:-Elements:-Shutdown())
                ),
                Maplets:-Elements:-BoxRow(
                    Maplets:-Elements:-BoxColumn(
                        Maplets:-Elements:-Label("EXPLANATION - Step-by-step student guide"),
                        Maplets:-Elements:-TextBox['ResultBox'](
                            'editable'=false, 'wrapped'=true, 'width'=46,'height'=34
                        )
                    ),
                    Maplets:-Elements:-BoxColumn(
                        Maplets:-Elements:-Label("PLOT: Z integer points / R real conic / C complex projection"),
                        Maplets:-Elements:-Plotter['Graph'](
                            'width'=740,'height'=530
                        )
                    )
                )
            )
        ),
        Maplets:-Elements:-Action['exampleHyperbolaAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_ExampleHyperbola")
        ),
        Maplets:-Elements:-Action['exampleEllipseAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_ExampleEllipse")
        ),
        Maplets:-Elements:-Action['exampleParabolaAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_ExampleParabola")
        ),
        Maplets:-Elements:-Action['examplePellAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_ExamplePell")
        ),
        Maplets:-Elements:-Action['exampleLinesAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_ExampleTwoLines")
        ),
        Maplets:-Elements:-Action['inputHelpAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_InputHelp")
        ),
        Maplets:-Elements:-Action['startupAction'](
            Maplets:-Elements:-RunWindow('MainWindow')
        ),
        Maplets:-Elements:-Action['bridgeAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_MathBridge")
        ),
        Maplets:-Elements:-Action['runAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_ShowReport")
        ),
        Maplets:-Elements:-Action['complexCurveAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_ComplexCurveReport"),
            Maplets:-Elements:-Evaluate('target'='Graph',
                                       'function'="GQS_M_ComplexCurvePlot")
        ),
        Maplets:-Elements:-Action['complexAction'](
            Maplets:-Elements:-Evaluate('function'="GQS_M_ComplexReport"),
            Maplets:-Elements:-Evaluate('target'='Graph',
                                          'function'="GQS_M_ComplexPlot")
        ),
        Maplets:-Elements:-Action['plotAction'](
            Maplets:-Elements:-Evaluate('target'='Graph',
                                          'function'="GQS_M_Plot")
        )
    );
    Maplets:-Display(M);
    return NULL;
end proc:

# Opens the interactive application after loading the file.
# Comment this out if you want to load the code without opening the window.
GQS_Maplet();

correction
 

Apparently, the LLM model is simply too small to be able to reason this out?

Download torus-clairaut-worksheet_29-9-2026.mw
torus-clairaut-worksheet-VERSION_B_29-9-2026.mw

I used the geom3d package 

restart:
with(geom3d):

# Twee vlakken
plane(P1, 2*x-y+z=3, [x,y,z]): plane(P2, x+2*y-z=-1, [x,y,z]):

# Punten op P1
point(A, 1, 0, 1): point(B, 0, 0, 3):

# Punten op P2
point(C, 0, 0, 1): point(E, 1, 0, 2):

draw([P1(color=blue,transparency=0.35,style=surface),
      P2(color=green,transparency=0.35,style=surface),
      A(color=red,symbol=solidsphere,symbolsize=18),
      B(color=red,symbol=solidsphere,symbolsize=18),
      C(color=red,symbol=solidsphere,symbolsize=18),
      E(color=red,symbol=solidsphere,symbolsize=18)],
      view=[-4..4,-4..4,-4..6],axes=boxed,labels=[x,y,z],
      scaling=constrained,orientation=[60,65]);

It's complicated to make a Geom3D package completely easy to use.
I've asked for outside help to point me in the right direction.

@Angelo Melino 

@Ronan 
A start module 
weinburger_-module_-maple_primesDEF_15-8-2026.mw

  The integral outcome

 for 0< a <Pi

restart:
with(codegen):

F := proc(x,y)
    x*y^2
end proc:

G_forward := GRADIENT(F, mode=forward):

G_forward;
                           G_forward


G_forward(2,3);
                             9, 12

foo := proc(params)
local b,c,d;

b := rhs(params[1]);
c := rhs(params[2]);
d := rhs(params[3]);

print(b,c,d);

end proc;

foo([b = 2, c = 3, d = 4]);
foo([b = 2, c = 3, d = 4]);
                            2, 3, 4

                            2, 3, 4

Once the procedure has finished, these local variables are no longer in use.
No global assignments were made.

Therefore, the second call is exactly the same as the first:

in this format : as fiirst order recurrence equation

 Alfred_F 590
After series option in dsolve  : order= 10 as example 

restart;

infolevel[dsolve] := 5:

ode := (1+f(x))*diff(f(x),x$2)=1+x:
ics := {f(0)=1,D(f)(0)=0}:

T := time():

sol := dsolve(ics union {ode}, f(x));

time()-T;

Does the AI in Maple itself provide better answers than when using AI outside of Maple?

I've been tinkering with the tensor a bit


 

experiment with different conditions for heated rod ( easy experiment, because there can be a lot more involved)   )

heat_pde_procedure_mprimes_16-7-2026.mw

use display command with ; 
The code was otherwise to large for uploading to the mprime server 
 

Download conducting_rod_in_magnetisch_veld_mprimesVersie_A_8-7-2026.mw