(*^ ::[paletteColors = 128; showRuler; currentKernel; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e8, 24, "Times"; ; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e6, 18, "Times"; ; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, L1, e6, 14, "Times"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, L1, a20, 18, "Times"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, L1, a15, 14, "Times"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, L1, a12, 12, "Times"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Times"; ; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 10, "Times"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, L1, 12, "Courier"; ; fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, L1, 10, "Geneva"; ; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Times"; ; fontset = Left Header, inactive, 12, "Times"; ; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, L1, 12, "Times"; ; fontset = Left Footer, inactive, 12, "Times"; ; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 10, "Times"; ; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Times"; ; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Times"; ; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Times"; ; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Times"; ; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Times"; ; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Times"; ; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Times"; ;] :[font = title; inactive; dontPreserveAspect; startGroup; ] Elliptic Curve Calculator :[font = text; inactive; dontPreserveAspect; center; ] Version 2.1 Ñ September 15th, 1992 Joseph H. Silverman Mathematics Department Brown University Providence, RI 02912 Email: ma420000@brownvm.bitnet Paul A. van Mulbregt Dragon Systems 320 Nevada Street Newton, MA 02160 Email: paulvm@dragonsys.com ;[s] 1:0,1;253,-1; 2:0,13,9,Times,0,12,0,0,0;1,16,12,Times,1,14,0,0,0; :[font = print; inactive; wordwrap; dontPreserveAspect; endGroup; ] A collection of Mathematica routines to perform calculations on elliptic curves defined over the rational numbers. Elliptic Curve Calculator copyright 1990-1992, Joseph H. Silverman. Mathematica is a trademark of Wolfram Research Inc. ;[s] 5:0,0;16,1;27,0;184,1;195,0;236,-1; 2:3,12,10,Courier,0,12,0,0,0;2,12,10,Courier,2,12,0,0,0; :[font = section; inactive; Cclosed; preserveAspect; startGroup; ] Implementation :[font = subsection; inactive; initialization; Cclosed; preserveAspect; startGroup; ] Setup the package context :[font = input; initialization; preserveAspect; endGroup; ] *) BeginPackage["EllipticCurves`EllipticCurveCalc`"] (* :[font = subsection; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Usage Messages :[font = text; inactive; Cclosed; preserveAspect; startGroup; ] Version and StartUp Message :[font = input; initialization; dontPreserveAspect; endGroup; ] *) ClearAll[$EllipticCurveCalcVersion, $EllipticCurveCalcVersionNumber]; $EllipticCurveCalcVersion="EllipticCurveCalculator v2.1"; $EllipticCurveCalcVersionNumber=2.1; EllipticCurveCalc::StartupNotice= "Elliptic Curve Calculator version `` \n Written by Joseph Silverman and Paul van Mulbregt.\n Copyright Joseph H. Silverman 1989-1992"; Message[EllipticCurveCalc::StartupNotice,$EllipticCurveCalcVersionNumber] (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Parameters and Miscellaneous Low-Level Functions :[font = input; initialization; dontPreserveAspect; endGroup; ] *) ClearAll[Origin, DiscriminantFactored, j, Discriminant,x,y,q]; Origin = P[Infinity,Infinity]; (* Zero point for curve *) DiscriminantFactored = Null; (* List of factors *) j:=j; Discriminant:=Discriminant; ord::usage="ord[n,p] computes the order of p dividing n."; RationalQ::usage="RationalQ[x] returns True if x is a rational number, False otherwise."; ClearECC::usage= "ClearECC[] clears some of the Elliptic Curve Calculator. Useful for an repetitive task involving many elliptic curves, which would tend to overburden the stack."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Entering and Displaying an Elliptic Curve :[font = input; initialization; dontPreserveAspect; endGroup; ] *) EllipticCurve::usage = "EllipticCurve[a1,a2,a3,a4,a6]: Enter generalized Weierstrass coefficients a1,a2,a3,a4,a6 \n Options: \n\t AlwaysFactorDiscriminant->True (default) or False. \n\t(Set to False to save time if calculations will not need the factorization of the discriminant. E.g. Graphing or adding points.) \n\t DisplayCurve->True (default) or False"; EllipticCurve::singular="\t Curve is singular !!"; EllipticCurve::notint= "The equation does not have integral coefficients."; Options[EllipticCurve] := {AlwaysFactorDiscriminant -> True, DisplayCurve -> True} EllipticCurveWNF::usage= "EllipticCurveWNF[a,b,c] enters the elliptic curve\n y^2=x^3+a*x^2+b*x+c"; OldCurve::usage = "OldCurve[n] reloads Curve n."; OldCurve::illnum="Curve number must be between 1 and ``."; DisplayEllipticCurve::usage = "DisplayEllipticCurve[] displays the Weierstrass equation in a nice format"; DisplayWeierstrassData::usage= "DisplayWeierstrassData[] displays a's, b's ,c's ,Discriminant, and j-invariant associated to the Weierstrass equation"; displayweierstrassdata::usage= "displayweierstrassdata[{a1,a2,a3,a4,a6}] displays a's, b's ,c's ,Discriminant, and j-invariant associated to the Weierstrass equation"; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Minimal Weierstrass Equations :[font = input; initialization; dontPreserveAspect; endGroup; ] *) MinimalQ::usage= "MinimalQ[] returns True if the Weierstrass equation is minimal, otherwise returns False. Does not affect the equation."; minimalQ::usage= "minimalQ[{a1,a2,a3,a4,a6}] returns True if the Weierstrass equation is minimal, otherwise returns False. Does not affect the equation."; MinimalEquation::usage= "MinimalEquation[] changes coordinates to produce a minimal Weierstrass equation in the standard form \n\t\t a[1] = 0,1, a[2] = -1,0,1, a[3] = 0,1. \n It returns a list {u,r,s,t} describing the change of coordinates used. To simultaneously change the coordinates of a list of points, use the command \n\t\t ShiftPoints[ {list of points}, MinimalEquation[] ]. \n To check if the equation is minimal without changing coordinates, use the command MinimalQ[]."; minimalequation::usage= "minimalequation[{a1,a2,a3,a4,a6}] returns a minimal Weierstrass equation in the standard form \n\t\t a[1] = 0,1, a[2] = -1,0,1, a[3] = 0,1. \n It returns a list {u,r,s,t} describing the change of coordinates used. \n To check if the equation is minimal without changing coordinates, use the command minimalQ[{a1,a2,a3,a4,a6}]."; ShiftPoints::usage= "ShiftPoints[ {P[x1,y1],P[x2,y2],...}, {u,r,s,t} ]: returns a list of the given points with their coordinates shifted by the indicated {u,r,s,t}."; LaskaAlgorithm::usage= "LaskaAlgorithm[] Returns a list {u,r,s,t} describing the change of coordinates needed to make the Weierstrass equation minimal and put it in the standard form \n\t\t a[1] = 0,1, a[2] = -1,0,1, a[3] = 0,1. \n In particular, u == 1 if and only if the equation is already minimal."; laskaalgorithm::usage= "laskallgorithm[{a1,a2,a3,a4,a6}] Returns a list {u,r,s,t} describing the change of coordinates needed to make the Weierstrass equation minimal and put it in the standard form \n\t\t a[1] = 0,1, a[2] = -1,0,1, a[3] = 0,1. \n In particular, u == 1 if and only if the equation is already minimal."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Adding, Subtracting, and Multiplying Points :[font = input; initialization; dontPreserveAspect; endGroup; ] *) ClearAll[P]; P::usage= "P[x,y] is a point on the curve with coordinates (x,y). The zero element may be denoted either by Origin or P[Infinity,Infinity]. Points may be added, subtracted, and multiplied by integers with the usual notation. P[x] will find a point on the curve with a given x-coordinate."; addpoints::usage= "addpoints[{a1,a2,a3,a4,a6}, P[x1,y1], P[x2,y2]] adds the two points on the curve." negativepoint::usage= "negativepoint[{a1,a2,a3,a4,a6}, P[x,y]] returns the negative of the point on the curve." PowerPoint::usage= "PowerPoint[{a1,a2,a3,a4,a6}, n, P[x,y]] returns n*P[x,y] computed on the curve." p::usage= "p[{a1,a2,a3,a4,a6},x] will find a point on the curve with a given x-coordinate."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Points of Finite Order :[font = input; initialization; dontPreserveAspect; endGroup; ] *) TorsionGroup::usage= "TorsionGroup[] returns the set of points of finite order on the curve."; TorsionGroup::usage= "torsiongroup[{a1,a2,a3,a4,a6}] returns the set of points of finite order on the curve."; PointOrder::usage= "PointOrder[P[x,y]]: returns the order of the given point, either as an integer or as Infinity."; pointorder::usage= "pointorder[{a1,a2,a3,a4,a6},P[x,y]]: returns the order of the given point, either as an integer or as Infinity."; MiscellaneousPoints::usage= "MiscellaneousPoints[] return a list of rational points found while computing other aspects of the curve, most notably (currently only!) points discovered while computing the TorsionGroupInfo[] for the curve." WhichTorsion::usage= "WhichTorsion[g] returns the type of the group of points g. g is optional. Runs faster when g has been computed." TorsionQ::usage= "TorsionQ[P[x,y]] returns True if P[x,y] is a torsion point, False otherwise." torsionQ::usage= "torsionQ[{a1,a2,a3,a4,a6}, P[x,y]] returns True if P[x,y] is a torsion point, False otherwise." TorsionGroup::usage= "TorsionGroup[] returns the set of points of finite order on the curve together with their orders, and the type of the abelian group." torsiongroup::usage= "torsiongroup[{a1,a2,a3,a4,a6}] returns the set of points of finite order on the curve together with their orders, and the type of the abelian group." TorsionGroupInfo::usage= "TorsionGroupInfo[{a1,a2,a3,a4,a6},] returns the set of points of finite order on the curve."; torsiongroupinfo::usage= "torsiongroupinfo[{a1,a2,a3,a4,a6}] returns the set of points of finite order on the curve."; TorsionGroupPoints::usage= "TorsionGroupPoints[] returns the list of rational Torsion points."; torsiongrouppoints::usage= "torsiongrouppoints[{a1,a2,a3,a4,a6}] returns the list of rational Torsion points."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Point Relations, Searching for Points :[font = input; initialization; dontPreserveAspect; endGroup; ] *) LatticeRelation::usage= "LatticeRelation[matrix] returns a relation using small integer coefficients approximately satisfied by the columns of the real valued matrix."; PointRelation::usage= "PointRelation[P[x1,x2],...] returns a linear relation satisfied by the given points. It returns the string \"Independent\" if the points are independent."; pointrelation::usage= "pointrelation[{a1, a2, a3, a4, a6},P[x1,x2],...] returns a linear relation satisfied by the given points. It returns the string \"Independent\" if the points are independent."; FindPoints::usage= "FindPoints[bound]: Search for rational points with height less than a given bound. \n If bound = {A,D}, search for points whose x-coordinate is a/d with |a| =< A and 1 =< d =< D. If bound = H is a number, same as FindPoints[{H,H}]."; findpoints::usage= "findpoints[{a1, a2, a3, a4, a6}bound]: Search for rational points with height less than a given bound. \n If bound = {A,D}, search for points whose x-coordinate is a/d with |a| =< A and 1 =< d =< D. If bound = H is a number, same as findpoints[{a1, a2, a3, a4, a6},{H,H}]."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Computing Local and Global Canonical Heights :[font = input; initialization; dontPreserveAspect; endGroup; ] *) LocalHeight::usage= "LocalHeight[P[x,y],p] computes the local height of the point P[x,y] for the place p. p = Infinity gives the archimedean contribution. For finite p, the Log[p] factor is not included, the contribution from the formal group is not included, and it is assumed that the Weierstrass equation is minimal at p."; localheight::usage= "localheight[{a1,a2,a3,a4,a6},P[x,y],p] computes the local height of the point P[x,y] for the place p. p = Infinity gives the archimedean contribution. For finite p, the Log[p] factor is not included, the contribution from the formal group is not included, and it is assumed that the Weierstrass equation is minimal at p."; Height::usage= "Height[P[x,y]]: The canonical height of P[x,y] with respect to the divisor (O)."; h::usage= "h[P[x,y]]: The canonical height of P[x,y] with respect to the divisor (O)."; h::nonminimal="Equation is not minimal."; height::usage= "height[{a1,a2,a3,a4,a6},P[x,y]]: The canonical height of P[x,y] with respect to the divisor (O)."; HeightPairing::usage= "HeightPairing[ P[x1,y1], P[x2,y2] ]: The canonical height pairing of the two points, normalized so that h[P] equals HeightPairing[P,P]."; heightpairing::usage= "heightpairing[{a1,a2,a3,a4,a6}, P[x1,y1], P[x2,y2] ]: The canonical height pairing of the two points, normalized so that height[{a1,a2,a3,a4,a6},P] equals heightpairing[{a1,a2,a3,a4,a6},P,P]."; HeightAngle::usage= "HeightAngle[ P[x1,y1], P[x2,y2] ]: The angle between the two points relative to the canonical height pairing."; heightangle::usage= "heightangle[{a1,a2,a3,a4,a6}, P[x1,y1], P[x2,y2] ]: The angle between the two points relative to the canonical height pairing."; HeightRegulator::usage= "HeightRegulator[P[x1,y1],P[x2,y2],...]: The height regulator for the given set of points. (N.B. This differs from the regulator in the Birch-Swinnerton-Dyer conjecture by a power of 2.)"; heightregulator::usage= "heightregulator[{a1,a2,a3,a4,a6},P[x1,y1],P[x2,y2],...]: The height regulator for the given set of points. (N.B. This differs from the regulator in the Birch-Swinnerton-Dyer conjecture by a power of 2.)"; HeightMatrix::usage= "HeightMatrix[P[x1,y1],P[x2,y2],...]: The matrix of height pairings for the given set of points."; heightmatrix::usage= "heightmatrix[{a1,a2,a3,a4,a6},P[x1,y1],P[x2,y2],...]: The matrix of height pairings for the given set of points."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Computing the Conductor and Reduction Types :[font = input; initialization; dontPreserveAspect; endGroup; ] *) ReductiontypeList::usage= "ReductiontypeList[prime] returns a list with the following information: \n\t\t the prime p \n\t\t order of the discriminant at p \n\t\t order of the conductor at p \n\t\t Kodaira reduction type (as a string) \n\t\t Birch--Swinnerton-Dyer \"fudge factor\" \n\t\t Good, Singular, Additive, Split or Non-Split, \n\t\t the coded reduction type."; reductiontypelist::usage= "reductiontypelist[{a1,a2,a3,a4,a6},prime] returns a list with the following information: \n\t\t the prime p \n\t\t order of the discriminant at p \n\t\t order of the conductor at p \n\t\t Kodaira reduction type (as a string) \n\t\t Birch--Swinnerton-Dyer \"fudge factor\" \n\t\t Good, Singular, Additive, Split or Non-Split, \n\t\t the coded reduction type."; Reductiontype::usage= "Reductiontype[prime] returns the reduction type of the elliptic curve in the following encoded form: \n 0 \t\t\t I(0) (good reduction) \n 2,3,4 \t\t II,III,IV \n -2,-3,-4 \t II*,III*,IV* \n 10 + n \t I(n), n>=1 \n -10 - n \t I(n)* \n 1 \t\t\t not minimal \n -1 \t\t singular"; reductiontype::usage= "reductiontype[{a1,a2,a3,a4,a6},prime] returns the reduction type of the elliptic curve in the following encoded form: \n 0 \t\t\t I(0) (good reduction) \n 2,3,4 \t\t II,III,IV \n -2,-3,-4 \t II*,III*,IV* \n 10 + n \t I(n), n>=1 \n -10 - n \t I(n)* \n 1 \t\t\t not minimal \n -1 \t\t singular"; ReductiontypeList::nonminimal= "The given curve is not minimal at the prime ``.\n\r Try with \n\t ``."; ConductorList::usage= "ConductorList[] displays a table giving primes of bad reduction, the exponent of the discriminant, the exponent of the conductor, and the Kodaira reduction type. To get the value of the conductor, use the command \"Conductor[]\"."; conductorlist::usage= "conductorlist[{a1,a2,a3,a4,a6},] displays a table giving primes of bad reduction, the exponent of the discriminant, the exponent of the conductor, and the Kodaira reduction type. To get the value of the conductor, use the command \"Conductor[]\"."; ConductorListFull::usage= "ConductorListFull[] displays a table giving primes of bad reduction, the exponent of the discriminant, the exponent of the conductor, the value of the \"fudge-factor\" c_p, and the Kodaira reduction type. To get the value of the conductor, use the command \"Conductor[]\"."; conductorlistfull::usage= "conductorlistfull[{a1,a2,a3,a4,a6}] displays a table giving primes of bad reduction, the exponent of the discriminant, the exponent of the conductor, the value of the \"fudge-factor\" c_p, and the Kodaira reduction type. To get the value of the conductor, use the command \"Conductor[]\"."; Conductor::usage= "Conductor[] returns the conductor of the elliptic curve."; conductor::usage= "conductor[{a1,a2,a3,a4,a6}] returns the conductor of the elliptic curve."; SingularPoint::usage= "SingularPoint[prime]: Returns a point P[x,y] which is a singular point on the curve modulo the indicated prime. If the curve has good reduction, the Origin is returned."; singularpoint::usage= "singularpoint[{a1,a2,a3,a4,a6},prime]: Returns a point P[x,y] which is a singular point on the curve modulo the indicated prime. If the curve has good reduction, the Origin is returned."; ChangeCoordinates::usage= "ChangeCoordinates[u,r,s,t]: Change Weierstrass coordinates by setting: x' = u^2 x + r, y' = u^3 y + s u^2 x + t, and recompute the a[i] coefficients. Does not recompute the other Weierstrass coefficients or the discriminant."; changecoordinates::usage= "changecoordinates[{a1,a2,a3,a4,a6},u,r,s,t]: Change Weierstrass coordinates by setting: x' = u^2 x + r, y' = u^3 y + s u^2 x + t. Returns the new a coefficients"; ChangeCoordinatesFull::usage= "ChangeCoordinatesFull[u,r,s,t] makes a full change of coordinates so that the a's, b's and c's are all changed."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Number Theory Functions :[font = input; initialization; preserveAspect; endGroup; ] *) SquareFreeQ::usage = "SquareFreeQ[n] is True if n is a square free integer, False otherwise."; FundamentalDiscriminantQ::usage = "FundamentalDiscriminantQ[n] is True if n is a fundamental discrimant of a quadratic field, False otherwise."; FundamentalDiscriminant::usage = "FundamentalDiscriminant[n] returns the fundamantal discriminant of the field Q(Sqrt[n])."; KroneckerSymbol::notcomp= "Unable to compute KroneckerSymbol as `` is not a valid discriminant."; KroneckerSymbol::usage= "KroneckerSymbol[d,n] returns the KroneckerSymbol of d, and n, where d is the discriminant of a quadratic field"; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Quadratic Twists :[font = input; initialization; preserveAspect; endGroup; ] *) QuadraticTwist::usage= "QuadraticTwist[d] returns a minimal equation for the quadratic twist of the elliptic curve by d.\n QuadraticTwist[d,n] returns a minimal equation for the quadratic twist of OldCurve[n] by d.\n QuadraticTwist[d,{a1,a2,a3,a4,a6}] returns a minimal equation for the quadratic twist of the curve {a1,a2,a3,a4,a6} by d.\n In all cases, d must be a fundamental quadratic discriminant." quadratictwist::usage= "quadratictwist[{a1,a2,a3,a4,a6},d] returns a minimal equation for the quadratic twist of the elliptic curve by d." QuadraticTwist::notfunddiscrimant= "The discrimant `` is not a fundamental quadratic discriminant."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Modp :[font = input; initialization; dontPreserveAspect; endGroup; ] *) FrobeniusTrace::usage= "FrobeniusTrace[prime] returns the trace of Frobenius for the given prime. If the curve has bad reduction, returns -1,1,0 if reduction is split, non-split, additive respectively. In all cases, this is the p th coefficient of the L-series."; frobeniustrace::usage= "frobeniustrace[{a1,a2,a3,a4,a6},prime] returns the trace of Frobenius for the given prime. If the curve has bad reduction, returns -1,1,0 if reduction is split, non-split, additive respectively. In all cases, this is the p th coefficient of the L-series."; SquareRootArray::usage="SquareRootArray[p] returns an array of length p-1, and the k-th entry is one of the two square roots of k, or -1 if no such square root exists mod p."; PointsModp::usage="PointsModp[p] returns a list of the solutions to the Elliptic Curve mod p. If you just want the number of solutions mod p, then the speedier function is\n p+1-LseriesCoefficient[p]."; pointsmodp::usage="pointsmodp[{a1,a2,a3,a4,a6},p] returns a list of the solutions to the Elliptic Curve mod p. If you just want the number of solutions mod p, then the speedier function is\n p+1-lseriescoefficient[{a1,a2,a3,a4,a6},p]."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Periods :[font = input; initialization; dontPreserveAspect; endGroup; ] *) RealPeriod::usage="RealPeriod[] computes the period that appears in the Birch-Swinnerton-Dyer conjecture. If the real locus is connected, this is the usual real period; if the real locus is disconnected, it is twice the usual real period."; realperiod::usage="realperiod[{a1,a2,a3,a4,a6}] computes the period that appears in the Birch-Swinnerton-Dyer conjecture. If the real locus is connected, this is the usual real period; if the real locus is disconnected, it is twice the usual real period."; ImaginaryPeriod::usage= "ImaginaryPeriod[] computes the imaginary period of the curve if the real locus is disconnected. It is the actual period, not twice it."; imaginaryperiod::usage= "imaginaryperiod[{a1,a2,a3,a4,a6}] computes the imaginary period of the curve if the real locus is disconnected. It is the actual period, not twice it."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] L-Series :[font = input; initialization; dontPreserveAspect; endGroup; ] *) FudgeFactor::usage= "FudgeFactor[] returns the Birch-Swinnerton-Dyer \"fudge factor\", equal to the product over the bad primes of the number of components on the Neron model."; fudgefactor::usage= "fudgefactor[{a1,a2,a3,a4,a6},] returns the Birch-Swinnerton-Dyer \"fudge factor\", equal to the product over the bad primes of the number of components on the Neron model."; LSeriesCoefficient::usage= "LSeriesCoefficient[n] returns the n th coefficient of the L-series of the curve." lseriescoefficient::usage= "lseriescoefficient[{a1,a2,a3,a4,a6},n] returns the n th coefficient of the L-series of the curve." LSeries::usage= "LSeries[s, options] returns the value of the L-series. \n\t Options: \n\t\t NumberOfTerms -> number of terms to use. (Default is 20.)" lseries::usage= "lseries[{a1,a2,a3,a4,a6},s, options] returns the value of the L-series. \n\t Options: \n\t\t NumberOfTerms -> number of terms to use. (Default is 20.)" Options[LSeries] = {NumberOfTerms -> 20}; FunctionalEquationSign::usage= "FunctionalEquationSign[n]: returns the sign of the functional equation of the L-series. For additive reduction uses a computation that takes n steps. If n is omitted, 20 steps are used." functionalequationsign::usage= "functionalequationsign[{a1,a2,a3,a4,a6},n]: returns the sign of the functional equation of the L-series. For additive reduction uses a computation that takes n steps. If n is omitted, 20 steps are used." FunctionalEquationSign::notcomp= "Unable to compute the sign of the functional equation. Try using FunctionalEquationSign[n] or LSeries[s,NumberOfTerms->n] with n larger than ``."; functionalequationsign::notcomp= "Unable to compute the sign of the functional equation. Try using functionalequationsign[{a1,a2,a3,a4,a6},n] or lseries[{a1,a2,a3,a4,a6},s,NumberOfTerms->n] with n larger than ``."; FunctionalEquationSignTwist::usage= "FunctionalEquationSignTwist[epsilon,Conductor,d,opts___] returns the sign of the functional equation of E(d) given the conductor and sign of the functional equation of E. \n If d is coprime to the conductor of E, then opts are ignored.\n If d is not coprime to the conductor of E, then the curve is twisted by d, the sign computed, then the old curve reloaded." functionalequationsigntwist::usage= "functionalequationsigntwist[{a1,a2,a3,a4,a6},d,opts___] returns the sign of the functional equation of E(d) given the conductor and sign of the functional equation of E. \n If d is coprime to the conductor of E, then opts are ignored.\n If d is not coprime to the conductor of E, then the curve is twisted by d, the sign computed" FunctionalEquationSignTwist::notfundisc= "Requires a fundamental discriminant." BSD::usage= "BSD[P[x1,y1],P[x2,y2]...] returns the conjectured leading term of the LSeries at s=1, if P[x1,y1],P[x2,y2]... is a basis for E(Q)/torsion, without the contribution from Shah.\n Conjecturally\n \t\t LSeries[s] \n \t\t limit ---------- = BSD[P[x1,y1],P[x2,y2],..P[xr,yr]] * Shah \n \t\t s->1 (s-1)^r\n and\n \t\tLSeries[1] = BSD[] * Shah if the rank is zero." bsd::usage= "bsd[{a1,a2,a3,a4,a6},P[x1,y1],P[x2,y2]...] returns the conjectured leading term of the LSeries at s=1, if P[x1,y1],P[x2,y2]... is a basis for E(Q)/torsion, without the contribution from Shah.\n Conjecturally\n \t\t lseries[{a1,a2,a3,a4,a6},s] \n \t\t limit --------------------------- = bsd[{a1,a2,a3,a4,a6},P[x1,y1],P[x2,y2],..P[xr,yr]] * Shah \n \t\t s->1 (s-1)^r\n and\n \t\t lseries[{a1,a2,a3,a4,a6},1] = bsd[{a1,a2,a3,a4,a6},] * Shah if the rank is zero." (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Modular Forms :[font = input; initialization; dontPreserveAspect; endGroup; ] *) PullBackLSeries::usage= "PullBackLSeries[n] returns the q-expansion of the modular form associated with the L-Series, up to powers of q^n."; pullbacklseries::usage= "pullbacklseries[{a1,a2,a3,a4,a6},n] returns the q-expansion of the modular form associated with the L-Series, up to powers of q^n."; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Graphing :[font = input; initialization; dontPreserveAspect; endGroup; ] *) PlotEllipticCurve::usage= "PlotEllipticCurve[xmax,options]: Plots the elliptic curve with x-values up to xmax. The usual options for Plot may be specified. If xmax is omitted, the default range is x < 5. The default label for the plot is the equation of the curve; to omit the label, specify the option PlotLabel -> None."; plotellipticcurve::usage= "plotellipticcurve[{a1,a2,a3,a4,a6},xmax,options]: Plots the elliptic curve with x-values up to xmax. The usual options for Plot may be specified. If xmax is omitted, the default range is x < 5. The default label for the plot is the equation of the curve; to omit the label, specify the option PlotLabel -> None."; Options[PlotEllipticCurve] := {PlotLabel -> DisplayEllipticCurve[]}; GraphPoints::usage= "GraphPoints[{P[x1,y1],P[x2,y2],...},options]: Create a graphics object displaying the indicated list of points. \n Options: \n\t PointLabels -> Automatic (default) or None." Options[GraphPoints] := {PointLabels->Automatic}; GraphTangents::usage= "GraphTangents[{P[x1,y1],P[x2,y2],...},options]: Create a graphics object displaying the tangent lines to the elliptic curve at the indicated points. \n Options: \n\t dx -> x-span of tangent line (default = 5) \n\t dy -> y-span of tangent line if vertical (default = 5) \n\t DisplayPoints -> True (default) or False \n\t PointLabels -> Automatic (default) or None"; Options[GraphTangents] := { dx -> 5, dy -> 5, DisplayPoints -> True }; graphtangents::usage= "graphtangents[{a1,a2,a3,a4,a6},{P[x1,y1],P[x2,y2],...},options]: Create a graphics object displaying the tangent lines to the elliptic curve at the indicated points. \n Options: \n\t dx -> x-span of tangent line (default = 5) \n\t dy -> y-span of tangent line if vertical (default = 5) \n\t DisplayPoints -> True (default) or False \n\t PointLabels -> Automatic (default) or None"; GraphSecants::usage= "GraphSecants[{P[x1,y1],P[x2,y2]},...},options]: Create a graphics object displaying the secant lines connecting the indicated pairs of points. \n Options: \n\t Extension -> 0.1 (default). Amount that lines extend beyond points (as a fraction of the length of the line.) \n\t DisplayPoints -> True (default) or False \n\t PointLabels -> Automatic (default) or None"; Options[GraphSecants] := { Extension -> 0.1, DisplayPoints -> True }; (* :[font = text; inactive; initialization; Cclosed; preserveAspect; startGroup; ] Equations of Secants and Tangents and the curve iteself. :[font = input; initialization; preserveAspect; endGroup; ] *) TangentLine::usage= "TangentLine[P[x,y]] returns the equation of the tangent line to the point P[x,y]."; tangentline::usage= "tangentline[{a1,a2,a3,a4,a6},P[x,y]] returns the equation of the tangent line to the point P[x,y]."; SecantLine::usage= "SecantLine[P[x,y],P[z,w]] returns the equation of the secant line through the points P[x,y], P[z,w]."; EquationCurve::usage= "EquationCurve[x,y] returns the equation of the elliptic curve in a form suitable for passing to Solve." (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Legendre Elliptic Curves y^2=x(x-1)(x-lambda) :[font = input; initialization; dontPreserveAspect; endGroup; ] *) OmegaOne::usage= "OmegaOne[lambda]: returns the first period, w1(lambda)."; OmegaOne::notcomp= "Unable to compute OmegaOne[``]. Need lambda close to 0 and 1."; OmegaOne::sing= "Unable to compute OmegaOne[lambda] with lambda=``. lambda=0 or lambda=1 give singular curves."; OmegaTwo::usage= "OmegaTwo[lambda]: returns the second period, w2(lambda)."; OmegaTwo::notcomp= "Unable to compute OmegaTwo[lambda] with lambda=``. Need lambda close to 0 and 1``."; OmegaTwo::sing= "Unable to compute OmegaTwo[lambda] with lambda=``. lambda=0 or lambda=1 give singular curves."; EllipticCurveLegendre::usage= "EllipticCurveLegendre[lambda] enters the elliptic curve y^2=x(x-1)(x-lambda)"; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Display the a, b, and c coefficients, the Discriminant, and the j-invariant :[font = input; initialization; dontPreserveAspect; endGroup; endGroup; ] *) Displayacoeffs::usage= "Displayacoeffs[] prints the a coefficients."; Displaybcoeffs::usage= "Displaybcoeffs[] prints the b coefficients."; Displayccoeffs::usage= "Displayccoeffs[] prints the c coefficients."; DisplayDiscj::usage= "DisplayDiscj[] prints the Discriminant and j-invariant."; (* :[font = subsection; inactive; initialization; Cclosed; preserveAspect; startGroup; ] Set the private context :[font = input; initialization; preserveAspect; endGroup; ] *) Begin["`Private`"] (* :[font = subsection; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Function Definitions for Elliptic Curve Calculator :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Parameters and Miscellaneous Low-Level Functions :[font = input; initialization; dontPreserveAspect; endGroup; ] *) curvenumber = 0; (* reference number for curve *) ClearAll[astack]; (* stack for a[i] coefs *) (* storecurve : Store the a[i] coefs in a stack *) storecurve := ( curvenumber++; astack[curvenumber]={a[1],a[2],a[3],a[4],a[6]}; ) (* ord[n,p] : Compute the order of p dividing n *) ord[n_Integer,p_Integer]:=Module[{order=0,m=n}, If[p<=1,Return[Infinity]]; If[n==0,Return[Infinity], (While[IntegerQ[m=m/p], (order++)]; Return[order];) ] ]; ord[n_Rational,p_Integer]:=ord[Numerator[n],p]-ord[Denominator[n],p]; SetAttributes[ord,Listable]; mod[m_Integer,n_Integer]:=Mod[m,n]; mod[a_,n_Integer]:=Mod[Numerator[a] * PowerMod[Denominator[a],-1,n],n]; SetAttributes[mod,Listable]; (* RationalQ[x]: test if x is in the rational numbers *) RationalQ[x_] := If[ Head[x]==Rational || Head[x]==Integer, True,False,False] (* displayfactorlist[list]: Displays a factored number in a nice format *) displayfactorlist[list_] := Apply[Times, list /. { {-1,1} -> -1, {p_Integer, e_Integer} -> HoldForm[p]^e }]; ClearECC[]:= ( curvenumber = 0; (* reference number for curve *) ClearAll[astack]; (* stack for a[i] coeffs *) miscellaneouspoints$0={}; functionalequationsign$0=0; ); miscellaneouspoints$0={}; functionalequationsign$0=0; a[1]=a[2]=a[3]=a[4]=a[6]=b[2]=b[4]=b[6]=b[8]=c[4]=c[6]=0; Discriminant=0; DiscriminantFactored=Null; j=Infinity; (* :[font = text; inactive; initialization; Cclosed; preserveAspect; startGroup; ] Coeffs from lists (not for export) :[font = input; initialization; preserveAspect; endGroup; ] *) bbfromaalist[{a1_,a2_,a3_,a4_,a6_}]:= {a1^2+4a2, 2a4+a1 a3, a3^2+4a6, a1^2 a6 + 4a2 a6 - a1 a3 a4 + a2 a3^2 - a4^2}; ccfrombblist[{b2_,b4_,b6_,b8_}]:= {b2^2 - 24b4, -b2^3 + 36 b2 b4 - 216 b6}; ccfromaalist[{a1_,a2_,a3_,a4_,a6_}]:= ccfrombblist[bbfromaalist[{a1,a2,a3,a4,a6}]]; discfromcclist[{c4_,c6_}]:=(c4^3-c6^2)/1728; discfrombblist[{b2_, b4_, b6_, b8_}]:= -b2^2 * b8 -8*b4^3 -27b6^2 +9b2*b4*b6; discfromaalist[{a1_,a2_,a3_,a4_,a6_}]:= discfrombblist[bbfromaalist[{a1,a2,a3,a4,a6}]]; discriminant[{a1_,a2_,a3_,a4_,a6_}]:= discfrombblist[bbfromaalist[{a1,a2,a3,a4,a6}]]; discriminantfactored[{a1_,a2_,a3_,a4_,a6_}]:= discriminantfactored[{a1,a2,a3,a4,a6}]= FactorInteger[discfromaalist[{a1,a2,a3,a4,a6}]]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Entering and Displaying an Elliptic Curve :[font = input; initialization; dontPreserveAspect; endGroup; ] *) EllipticCurve[a1_,a2_,a3_,a4_,a6_,opts___] := ( Unprotect[a,b,c,j,Discriminant, DiscriminantFactored]; ClearAll[a,b,c]; a[1]=a1; a[2]=a2; a[3]=a3; a[4]=a4; a[6]=a6; storecurve; b[2]=a1^2+4*a2; b[4]=a1*a3+2*a4; b[6]=a3^2+4*a6; b[8]=a1^2*a6+4*a2*a6-a1*a3*a4+a2*a3^2-a4^2; c[4]=b[2]^2-24*b[4]; c[6]=-b[2]^3+36*b[2]*b[4]-216*b[6]; Discriminant=(c[4]^3-c[6]^2)/1728; miscellaneouspoints$0={}; functionalequationsign$0=0; If[AlwaysFactorDiscriminant/.{opts} /.Options[EllipticCurve], DiscriminantFactored=FactorInteger[Discriminant], DiscriminantFactored=Null, DiscriminantFactored=Null]; discriminantfactored[{a1,a2,a3,a4,a6}]=DiscriminantFactored; If[Discriminant==0, Message[EllipticCurve::singular], j=c[4]^3/Discriminant] ; Protect[a,b,c,j,Discriminant,DiscriminantFactored]; If[DisplayCurve/.{opts} /.Options[EllipticCurve], Print[StringForm["Curve ``: ",curvenumber], DisplayEllipticCurve[]] ]; ) EllipticCurve[{a1_,a2_,a3_,a4_,a6_},opts___]:= EllipticCurve[a1,a2,a3,a4,a6,opts]; EllipticCurveWNF[a_,b_,c_,opts___]:=EllipticCurve[0,a,0,b,c,opts]; OldCurve[n_] := If[ 1 <= n <= curvenumber, Apply[EllipticCurve,astack[n]], Message[OldCurve::illnum,curvenumber] ]; OldCurve[n_,opts__] := If[ 1 <= n <= curvenumber, Apply[EllipticCurve,Append[astack[n],opts]], Message[OldCurve::illnum,curvenumber] ]; DisplayEllipticCurve[] := displayellipticcurve[{a[1],a[2],a[3],a[4],a[6]}]; pm[n_,s_] := Which[ n==0,"", n==1,StringForm[" + ``",s], n==-1,StringForm[" - ``",s], n<0,StringForm[" - `` ``",-n,s], True,StringForm[" + `` ``",n,s] ]; displayellipticcurve[{a1_,a2_,a3_,a4_,a6_}]:= StringForm["`````` = ````````", "y^2",pm[a1,"xy"],pm[a3,"y"], "x^3",pm[a2,"x^2"],pm[a4,"x"], pm[Sign[a6],Abs[a6]] ] ; DisplayWeierstrassData[] := Block[{}, Print["a coefficients = ",{a[1],a[2],a[3],a[4],a[6]}]; Print[]; Print["b coefficients = ",{b[2],b[4],b[6],b[8]}]; Print[]; Print["c coefficients = ",{c[4],c[6]}]; Print[]; Print["Discriminant = ",Discriminant]; If[Length[DiscriminantFactored]>0, Print[" = ", displayfactorlist[DiscriminantFactored]]]; Print[]; Print["j invariant = ",j]; If[Length[DiscriminantFactored]>0, Print[" = ", displayfactorlist[FactorInteger[c[4]]]^3/ displayfactorlist[DiscriminantFactored]]]; ]; displayweierstrassdata[{a1_,a2_,a3_,a4_,a6},factordisc_:TRUE] := Module[{b2,b4,b6,b8,c4,c6}, Print["a coefficients = ",{a1,a2,a3,a4,a6}]; Print[]; {b2,b4,b6,b8}=bbfromaalist[{a1,a2,a3,a4,a6}]; Print["b coefficients = ",{b2,b4,b6,b8}]; Print[]; {c4,c6} = ccfrombblist[{b2,b4,b6,b8}]; Print["c coefficients = ",{c4,c6}]; Print[]; Print["Discriminant = ",discriminantfrombblist[{b2,b4,b6,b8}]]; If[factordisc, Print[" = ", displayfactorlist[discriminantfactored[{a1,a2,a3,a4,a6}]]]]; Print[]; Print["j invariant = ",j]; If[factordisc, Print[" = ", displayfactorlist[FactorInteger[c4]]^3/ displayfactorlist[discriminantfactored[{a1,a2,a3,a4,a6}]] ] ]; ]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Minimal Weierstrass Equations :[font = input; initialization; dontPreserveAspect; endGroup; ] *) MinimalQ[] := minimalQ[{a[1], a[2], a[3], a[4], a[6]}]; minimalQ[{a1_, a2_, a3_, a4_, a6_}]:= minimalQ[{a1, a2, a3, a4, a6}]= If[laskaalgorithm[{a1,a2,a3,a4,a6}][[1]]==1,True,False,False] minimalQ[a1_, a2_, a3_, a4_, a6_]:= minimalQ[{a1, a2, a3, a4, a6}]; minimalQ[xxx___]:=False; LaskaAlgorithm[] := laskaalgorithm[{a[1],a[2],a[3],a[4],a[6]}]; laskaalgorithm[{a1_,a2_,a3_,a4_,a6_}] := Module[ {u,r,s,t,umax,c4,c6}, (* Check equation is integral *) If[Max[Map[Denominator, {a1,a2,a3,a4,a6}]]>1, Message[EllipticCurve::notint]; Return[Null]]; {c4,c6}=ccfromaalist[{a1,a2,a3,a4,a6}]; (* Main routine *) umax=GCD[c4,c6]; If[Abs[umax]!=1, umax = Apply[Times, Map[ #^Floor[ Min[ord[c4,#]/4, ord[c6,#]/6]]&, Thread[FactorInteger[umax]][[1]] ] ] ]; u = 2^ord[umax,2]; umax = umax/u; While[ u>1 && Not[laskatest[u,{a1,a2,a3,a4,a6,c4,c6}][[1]]], u = u/2 ]; u = u*umax; While[ u>1 && Not[laskatest[u,{a1,a2,a3,a4,a6,c4,c6}][[1]]], u = u/3 ]; Return[Rest[laskatest[u,{a1,a2,a3,a4,a6,c4,c6}]]]; ] (* This is the test used in the middle of the Laska routine to find a minimal equation. It returns {False} if u is not a possible scaling factor, and {True,u,r,s,t} if {u,r,s,t} is a possible change of variables. *) laskatest[u_,{aa1_,aa2_,aa3_,aa4_,aa6_,c4_,c6_}] := Module[{a1,a2,a3,a4,a6,x,y,r,s,t}, x = c4/u^4; y = c6/u^6; a1 = Mod[x,2]; a2 = Switch[Mod[-a1-y,3],0,0,1,1,2,-1]; If[(a1==0 && !IntegerQ[y/8]) || (a1==1 && !IntegerQ[(x-1)/8]),Return[{False}]]; a3 = Which[ a1==0, Mod[y/8,2], a1==1, Mod[(x-1)/8+a2,2] ]; a4 = (-x + (a1+4a2)^2)/48 - a1*a3/2; a6 =(-y-(a1+4a2)^3+36(a1+4a2)(a1*a3+2a4))/864-a3/4; If[!IntegerQ[a4] || !IntegerQ[a6],Return[{False}]]; s = (u*a1-aa1)/2; r = (u^2*a2-aa2+s*aa1+s^2)/3; t = (u^3*a3-aa3-r*aa1)/2; If[!IntegerQ[s] || !IntegerQ[r] || !IntegerQ[t], Return[{False}] ]; If[u^4*a4 != aa4-s*aa3+2r*aa2-(t+r*s)aa1+3r^2-2s*t, Return[{False}] ]; If[u^6*a6 != aa6+r*aa4+r^2*aa2+r^3-t*aa3-t^2-r*t*aa1, Return[{False}] ]; Return[{True,u,r,s,t}] ]; MinimalEquation[] := Module[{res}, res=minimalequationandurst[{a1,a2,a3,a4,a6}]; EllipticCurve[Apply[Sequence,First[res]] ]; res[[2]] ] minimalequation[{a1_,a2_,a3_,a4_,a6_}] := minimalequationandurst[{a1,a2,a3,a4,a6}]; minimalequationandurst[{a1_,a2_,a3_,a4_,a6_}] := Module[ {u,r,s,t}, {u,r,s,t} = laskaalgorithm[{a1,a2,a3,a4,a6}]; {{(a1+2s)/u, (a2-s*a1+3r-s^2)/u^2, (a3+r*a1+2t)/u^3, (a4-s*a3+2r*a2-(t+r*s)a1+3r^2-2s*t)/u^4, (a6+r*a4+r^2*a2+r^3-t*a3-t^2-r*t*a1)/u^6}, {u,r,s,t}} ]; ShiftPoints[list:{P[_,_]..}, {u_,r_,s_,t_}] := list/. P[x_,y_] -> P[(x-r)/u^2,(y-s*(x-r)-t)/u^3]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Adding, Subtracting, and Multiplying Points :[font = input; initialization; dontPreserveAspect; endGroup; ] *) (* Set up addition/multiplication rules for points *) Format[P[Infinity,Infinity]]:=HoldForm[Origin]; (* Rules to Add Points *) P/: P[x1_,y1_] + P[x2_,y2_] := Module[ {Slope,x3}, If[ x1!=x2, Slope = (y2-y1)/(x2-x1), Slope = (3x1^2+2*a[2]*x1+a[4]-a[1]*y1)/ (2y1+a[1]*x1+a[3]), Slope = (y2-y1)/(x2-x1)]; P[x3=Slope^2+a[1]*Slope-a[2]-x1-x2, -Slope*(x3-x1)-y1-a[1]x3-a[3]] ] (* Special rule for adding origin *) P/: P[x_,y_] + Origin := P[x,y] (* Special rule for adding negatives *) P/: P[x1_,y1_] + P[x2_,y2_] := Origin /; x1==x2 && (y1!=y2 || 2y1+a[1]*x1+a[3]==0) (* Rules for Negating and Subtracting points *) P/: -P[x_,y_] := P[x,-(y+a[1]*x+a[3])] P/: -Origin := Origin P/: P[x1_,y1_] - P[x2_,y2_] := P[x1,y1] + (-P[x2,y2]) (* Rules for Multiplying points by integers *) (* n*P computed by recursive binary expansion of |n| *) P/: 0*P[x_,y_] := Origin P/: 1*P[x_,y_] := P[x,y] P/: 2*P[x_,y_] := P[x,y]+P[x,y] P/: n_Integer*P[x_,y_] := PowerPoint[n,P[x,y]] P/: PowerPoint[0,P[x_,y_]] := Origin P/: PowerPoint[1,P[x_,y_]] := P[x,y] P/: PowerPoint[n_,P[x_,y_]] := 2*PowerPoint[n/2,P[x,y]] /; n>1 && EvenQ[n] P/: PowerPoint[n_,P[x_,y_]] := 2*PowerPoint[(n-1)/2,P[x,y]] + P[x,y] \ /; n>1 && OddQ[n] P/: PowerPoint[n_,P[x_,y_]] := - PowerPoint[-n,P[x,y]] /; n<0 (* Find a point on the curve with a given x-coordinate *) P[x_] := P[x,y]/.Solve[y^2+a[1]x y+a[3]y -x^3-a[2]x^2-a[4]x-a[6]==0,y][[1]]; AddPoints[P[x1_,y1_], P[x2_,y2_]] := P[x1,y1] + P[x2,y2]; NegativePoint[P[x_,y_]]:=-P[x,y]; addpoints[P[x1_,y1_], P[x2_,y2_]] := P[x1,y1] + P[x2,y2]; negativepoint[P[x_,y_]]:=-P[x,y]; addpoints[{a1_, a2_, a3_, a4_, a6_}, P[x_,y_] , Origin] := P[x,y]; addpoints[{a1_, a2_, a3_, a4_, a6_},P[x1_,y1_],P[x2_,y2_]]:= Origin /; (x1==x2 && (y1!=y2 || 2y1+a1*x1+a3==0)) addpoints[{a1_, a2_, a3_, a4_, a6_},P[x1_,y1_],P[x2_,y2_]]:= Module[ {Slope,x3}, If[ x1!=x2, Slope = (y2-y1)/(x2-x1), Slope = (3x1^2+2*a2*x1+a4-a1*y1)/ (2y1+a1*x1+a3), Slope = (y2-y1)/(x2-x1)]; P[x3=Slope^2+a1*Slope-a2-x1-x2, -Slope*(x3-x1)-y1-a1*x3-a3] ]; negativepoint[{a1_, a2_, a3_, a4_, a6_},Origin]=Origin; negativepoint[{a1_, a2_, a3_, a4_, a6_},P[x_,y_]] := P[x,-(y+a1*x+a3)]; powerpoint[{a1_, a2_, a3_, a4_, a6_},0,P[x_,y_]] := Origin; powerpoint[{a1_, a2_, a3_, a4_, a6_},1,P[x_,y_]] := P[x,y]; powerpoint[{a1_, a2_, a3_, a4_, a6_},2,P[x_,y_]] := addpoints[{a1, a2, a3, a4, a6},P[x,y],P[x,y]]; powerpoint[{a1_, a2_, a3_, a4_, a6_},n_,P[x_,y_]] := powerpoint[{a1, a2, a3, a4, a6},2, powerpoint[{a1, a2, a3, a4, a6}, n/2, P[x,y]] ]/; n>1 && EvenQ[n]; powerpoint[{a1_, a2_, a3_, a4_, a6_},n_,P[x_,y_]] := addpoints[{a1, a2, a3, a4, a6}, powerpoint[{a1, a2, a3, a4, a6},2,powerpoint[{a1, a2, a3, a4, a6},(n-1)/2,P[x,y]]], P[x,y]] /; n>1 && OddQ[n]; powerpoint[{a1_, a2_, a3_, a4_, a6_},n_,P[x_,y_]] := negativepoint[{a1, a2, a3, a4, a6}, powerpoint[-n, P[x,y]]] /; n<0 powerpoint[a___] := powerpoint[a]; p[{a1_, a2_, a3_, a4_, a6_},x_] := P[x,y]/.Solve[y^2+a1*x y+a3*y-x^3-a2*x^2-a4*x-a6==0,y][[1]]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Points of Finite Order :[font = input; initialization; dontPreserveAspect; endGroup; ] *) PointOrder[P[x_,y_]] := Module[ {n,nP=P[x,y]}, If[ nP == Origin, Return[1], ]; Do[nP = nP+P[x,y]; If[ nP == Origin, Return[n], ]; If[ n==12, Return[Infinity]]; If[ Denominator[nP[[1]]]>4, Return[Infinity], ], {n,2,12} ] ]/; (RationalQ[x] && RationalQ[y]) PointOrder[P[x_,y_]]:=Infinity/;!(RationalQ[x] && RationalQ[y]) pointorder[{a1_, a2_, a3_, a4_, a6_},P[x_,y_]] := Module[ {n,nP=P[x,y]}, If[ nP == Origin, Return[1], ]; Do[nP = addpoints[{a1, a2, a3, a4, a6}, nP, P[x,y]]; If[ nP == Origin, Return[n], ]; If[ n==12, Return[Infinity]]; If[ Denominator[nP[[1]]]>4, Return[Infinity], ], {n,2,12} ] ]/; (RationalQ[x] && RationalQ[y]) pointorder[{a1_, a2_, a3_, a4_, a6_},P[x_,y_]] := Infinity/;!(RationalQ[x] && RationalQ[y]); whichgroup[tt_List]:= Switch[Sort[tt], {1},"Z1", {1,2},"Z2", {1,3,3},"Z3", {1,2,4,4},"Z4", {1,5,5,5,5},"Z5", {1,2,3,3,6,6},"Z6", {1,7,7,7,7,7,7},"Z7", {1,2,4,4,8,8,8,8},"Z8", {1,3,3,9,9,9,9,9,9},"Z9", {1,2,5,5,5,5,10,10,10,10},"Z10", {1,2,3,3,4,4,6,6,12,12,12,12},"Z12", {1,2,2,2}, "Z2 * Z2", {1,2,2,2,4,4,4,4}, "Z2 * Z4", {1,2,2,2,3,3,6,6,6,6,6,6}, "Z2 * Z6", {1,2,2,2,4,4,4,4,8,8,8,8,8,8,8,8},"Z2 * Z8", _,"Unknown" ] WhichTorsion[t_:{}]:= If[t=={}, TorsionGroupInfo[][[2]], whichgroup[Map[PointOrder,t]] ] TorsionQ[P[x_,y_]]:= Module[{eightP,twoP=2*P[x,y]}, If[Denominator[twoP[[1]]]>4,Return[False],]; eightP=4*twoP; If[eightP==Origin,Return[True],]; If[Denominator[eightP[[1]]]>4,Return[False],]; eightP=eightP+(8*eightP); If[(eightP==Origin)||(eightP==twoP), Return[True], Return[False], Return[False] ] ]/; (RationalQ[x] && RationalQ[y]) torsionQ[{a1_, a2_, a3_, a4_, a6_},P[x_,y_]]:= Module[{eightP,twoP=addpoints[{a1, a2, a3, a4, a6}, P[x,y], P[x,y]]}, If[Denominator[twoP[[1]]]>4,Return[False],]; eightP=powerpoint[{a1, a2, a3, a4, a6},4,twoP]; If[eightP==Origin, Return[True],]; If[Denominator[eightP[[1]]]>4, Return[False],]; eightP=powerpoint[{a1, a2, a3, a4, a6},9,twoP]; If[(eightP==Origin)||(eightP==twoP), Return[True], Return[False], Return[False] ] ]/; (RationalQ[x] && RationalQ[y]) TorsionQ[P[x_,y_]]:=False; TorsionGroupInfo[]:=torsiongroupinfowlist[{a[1],a[2],a[3],a[4],a[6]}]; torsiongroupinfo[{a1_,a2_,a3_,a4_,a6_}]:=torsiongroupinfowlist[{a1,a2,a3,a4,a6}]; torsiongroupinfowlist[{a1_, a2_, a3_, a4_, a6_}] := torsiongroupinfowlist[{a1, a2, a3, a4, a6}] = (* Save value *) Module[{listpts,miscpts,b2,b4,b6,b8}, {b2,b4,b6,b8}=bbfromaalist[{a1, a2, a3, a4, a6}]; listpts=Map[{#,pointorder[{a1,a2,a3,a4,a6},#]}&, Cases[ Flatten[ Map[ Thread[{x/.Solve[4x^3+b2*x^2+2*b4*x+b6==#^2,x],{#,#,#}}]&, Join[{0}, Divisors[ Apply[Times, discriminantfactored[{a1, a2, a3, a4, a6}]/. {p_Integer,e_Integer}->p^Floor[e/2] ] ] ] ], 1], {x_?RationalQ,y_} -> P[x,(y-a1*x-a3)/2] ] ]; listpts=Union[{{Origin,1}},listpts,listpts/.{exp1_,exp2_Integer}->{negativepoint[{a1,a2,a3,a4,a6},exp1],exp2}]; miscpts=Select[listpts,(!IntegerQ[#[[2]]])&]; If[(miscpts=={})||({a1,a2,a3,a4,a6}!={a[1],a[2],a[3],a[4],a[6]}),, miscellaneouspoints$0=Union[miscellaneouspoints$0,Thread[miscpts][[1]]] ]; listpts=Select[listpts,(IntegerQ[#[[2]]])&]; listpts=Append[{listpts},whichgroup[Thread[listpts][[2]]]] ] TorsionGroupPoints[]:=Thread[TorsionGroupInfo[][[1]]][[1]]; TorsionGroup[]:=TorsionGroupPoints[]; torsiongrouppoints[{a1_,a2_,a3_,a4_,a6_}]:= Thread[torsiongroupinfowlist[{a1,a2,a3,a4,a6}][[1]]][[1]]; torsiongroup[{a1_,a2_,a3_,a4_,a6_}]:=torsiongrouppoints[{a1,a2,a3,a4,a6}]; MiscellaneousPoints[]:= miscellaneouspoints$0; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Point Relations, Searching for Points :[font = input; initialization; dontPreserveAspect; endGroup; ] *) LatticeRelation[matrix_?MatrixQ]:= Drop[ LatticeReduce[ Thread[Join[ IdentityMatrix[Length[matrix]], Map[Round[#*10^6]&,matrix] ] ] ] [[1]], -Length[matrix] ] PointRelation[listofpoints__] := Module[ {htmat, relation}, htmat = HeightMatrix[listofpoints][[1]]; If[ Abs[Det[htmat]] > 10^-4, Return["Independent"] ]; relation = LatticeRelation[htmat]; If[ relation.{listofpoints} == Origin, Return[Inner[HoldForm[#1]*#2&, relation,{listofpoints}]==0 ] ]; Return["Could not determine dependence"] ] PointRelationTorsion[listofpoints__] := Module[{htmat, relation, pt}, htmat = HeightMatrix[listofpoints][[1]]; If[ Abs[Det[htmat]] > 10^-4, Return["Independent"] ]; relation = LatticeRelation[htmat]; pt = relation.{listofpoints}; If[ MemberQ[TorsionGroupPoints[],pt], Return[Inner[HoldForm[#1]*#2&, relation,{listofpoints}]==pt ] ]; Return["Could not determine dependence"] ] pointrelation[{a1_, a2_, a3_, a4_, a6_},listofpoints__] := Module[{htmat, relation}, htmat = heightmatrix[{a1, a2, a3, a4, a6},listofpoints][[1]]; If[ Abs[Det[htmat]] > 10^-4, Return["Independent"] ]; relation = LatticeRelation[htmat]; If[ relation.{listofpoints} == Origin, Return[Inner[HoldForm[#1]*#2&, relation,{listofpoints}]==0 ] ]; Return["Could not determine dependence"] ] pointrelationtorsion[{a1_, a2_, a3_, a4_, a6_},listofpoints__] := Module[{htmat, relation,pt}, htmat = heightmatrix[{a1, a2, a3, a4, a6},listofpoints][[1]]; If[ Abs[Det[htmat]] > 10^-4, Return["Independent"] ]; relation = LatticeRelation[htmat]; pt = relation.{listofpoints}; If[ MemberQ[TorsionGroupPoints[],pt], Return[Inner[HoldForm[#1]*#2&, relation,{listofpoints}]==pt ] ]; Return["Could not determine dependence"] ] FindPoints[bound_] := Module[{NumBound,DenBound,xmin,x,d,n,pointlist={}}, {NumBound,DenBound} = If[NumberQ[bound], {bound,Floor[N[Sqrt[bound]]]}, {bound[[1]],Floor[N[Sqrt[bound[[2]]]]]} ]; xmin = Min[Select[ (x/.{ToRules[ NRoots[4x^3+b[2]x^2+2b[4]x+b[6]==0,x]]}), Im[#]==0&]]; Do[ If[IntegerQ[Sqrt[4*n^3+b[2]*n^2*d^2+2*b[4]*n*d^4+b[6]*d^6]], AppendTo[pointlist,P[n/d^2]], ], {d,1,DenBound}, {n,Max[-NumBound,Floor[xmin*d^2]],NumBound} ]; Return[Union[pointlist,Minus[pointlist]]] ]; findpoints[{a1_, a2_, a3_, a4_, a6_},bound_] := Module[{NumBound,DenBound,xmin,x,d,n,pointlist={},b2,b4,b6,b8}, {NumBound,DenBound} = If[NumberQ[bound], {bound,Floor[N[Sqrt[bound]]]}, {bound[[1]],Floor[N[Sqrt[bound[[2]]]]]} ]; {b2,b4,b6,b8}=bbfromaalist[{a1, a2, a3, a4, a6}]; xmin = Min[Select[ (x/.{ToRules[ NRoots[4x^3+b2*x^2+2b4*x+b6==0,x]]}), Im[#]==0&]]; Do[ If[IntegerQ[Sqrt[4*n^3+b2*n^2*d^2+2*b4*n*d^4+b6*d^6]], AppendTo[pointlist,p[{a1, a2, a3, a4, a6},n/d^2]], ], {d,1,DenBound}, {n,Max[-NumBound,Floor[xmin*d^2]],NumBound} ]; Return[Union[pointlist,Minus[pointlist]]] ]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Computing Local and Global Canonical Heights :[font = input; initialization; dontPreserveAspect; endGroup; ] *) LocalHeight[P[x_,y_],p_Integer] := localheight[{a[1],a[2],a[3],a[4],a[6]}, P[x,y], p]; localheight[{a1_, a2_, a3_, a4_, a6_}, P[x_, y_], p_Integer]:= Module[ {c4,c6,b2,b4,b6,b8,AA,BB,CC,DD}, {b2,b4,b6,b8}=bbfromaalist[{a1, a2, a3, a4, a6}]; {c4,c6}=ccfrombblist[{b2,b4,b6,b8}]; AA=ord[Numerator[3x^2+2*a2*x+a4-a1*y],p]; BB=ord[Numerator[2y+a1*x+a3],p]; CC=ord[Numerator[3x^4+b2*x^3+3*b4*x^2+3b6*x+b8],p]; DD=ord[discfrombblist[{b2,b4,b6,b8}],p]; Return[ Which[ AA==0 || BB==0, 0, ord[c4,p]==0, Min[BB,DD/2]*(Min[BB,DD/2]-DD)/(2DD), CC>=3BB, -BB/3, True, -CC/8 ] ] ] LocalHeight[P[x_,y_],Infinity] := localheight[{a[1],a[2],a[3],a[4],a[6]}, P[x,y], Infinity] localheight[{a1_, a2_, a3_, a4_, a6_}, P[x_, y_], Infinity]:= Module[{bb,bp,w,z,t,OrigEq,lambda,mu,n}, bb=bbfromaalist[{a1, a2, a3, a4, a6}]; bp[2]=bb[[2/2]]-12; bp[4]=bb[[4/2]]-bb[[2/2]]+6; bp[6]=bb[[6/2]]-2bb[[4/2]]+bb[[2/2]]-4; bp[8]=bb[[8/2]]-3bb[[6/2]]+3bb[[4/2]]-bb[[2/2]]+3; If[Abs[x]>=.5, (t=1./x; OrigEq=True;), (t=1./(x+1.); OrigEq=False;) ]; lambda=-.5*Log[Abs[t]]; n=0; mu=0.; Do[ If[OrigEq, w=4.t+bb[[2/2]]t^2+2bb[[4/2]]t^3+bb[[6/2]]t^4; z=1.-bb[[4/2]]t^2-2.bb[[6/2]]t^3-bb[[8/2]]t^4; If[ Abs[w] <= 2 Abs[z], mu=mu+.25^n Log[Abs[z]]; t=w/z, mu=mu+.25^n Log[Abs[z+w]];t=w/(z+w); OrigEq=False ]; , w=4.t+bp[2]t^2+2bp[4]t^3+bp[6]t^4; z=1.-bp[4]t^2-2.bp[6]t^3-bp[8]t^4; If[ Abs[w] <= 2 Abs[z], mu=mu+.25^n Log[Abs[z]]; t=w/z, mu=mu+.25^n Log[Abs[z-w]];t=w/(z-w); OrigEq=True ]; ]; n=n+1, {Ceiling[ 5 Precision[1.]/3 +.5 + .75 Log[7.+(4/3)*Log[Max[Abs[{a1, a2, a3, a4, a6}]]]]] } ]; Return[N[lambda+.125*mu]]; ] h[Origin]:=0; h[P[x_,y_]] := hwlist[{a[1],a[2],a[3],a[4],a[6]},P[x,y]]; SetAttributes[h,Listable]; Height[P[x_,y_]] := hwlist[{a[1],a[2],a[3],a[4],a[6]},P[x,y]]; SetAttributes[h,Listable]; height[{a1_,a2_,a3_,a4_,a6_},P[x_,y_]] := hwlist[{a1,a2,a3,a4,a6},P[x,y]]; hwlist[{a1_,a2_,a3_,a4_,a6_},Origin]:=0; hwlist[{a1_,a2_,a3_,a4_,a6_},P[x_,y_]] := Module[{temp}, temp = Chop[ localheight[{a1, a2, a3, a4, a6},P[x,y],Infinity] + .5*N[Log[Denominator[x]]] + Map[localheight[{a1, a2, a3, a4, a6},P[x,y],#]&, Select[Thread[discriminantfactored[{a1, a2, a3, a4, a6}]][[1]], #>1& ] ] \ .N[Log[Select[Thread[discriminantfactored[{a1, a2, a3, a4, a6}]][[1]], #>1& ] ] ] ]; (*Save height for future reference*) Unprotect[hwlist]; hwlist[{a1, a2, a3, a4, a6},P[x,y]] = temp; Protect[hwlist]; Return[temp]; ]/;minimalQ[{a1, a2, a3, a4, a6}] hwlist[list_List,P[x,y]]:=Message[h::nonminimal]; HeightPairing[P[x1_,y1_],P[x2_,y2_]] := .5*(h[P[x1,y1]+P[x2,y2]]-h[P[x1,y1]]-h[P[x2,y2]]) heightpairing[{a1_,a2_,a3_,a4_,a6_},P[x1_,y1_],P[x2_,y2_]] := .5*(hwlist[{a1, a2, a3, a4, a6},addpoints[{a1, a2, a3, a4, a6},P[x1,y1], P[x2,y2]]] -hwlist[{a1, a2, a3, a4, a6},P[x1,y1]]- hwlist[{a1, a2, a3, a4, a6},P[x2,y2]]); HeightAngle[P[x1_,y1_],P[x2_,y2_]] := If[ h[P[x1,y1]]*h[P[x2,y2]]==0., 0., ArcCos[HeightPairing[P[x1,y1],P[x2,y2]]/ (Sqrt[h[P[x1,y1]]]*Sqrt[h[P[x2,y2]]]) ] ]; heightangle[{a1_,a2_,a3_,a4,_a6_},P[x1_,y1_],P[x2_,y2_]] := If[ hwlist[{a1, a2, a3, a4, a6},P[x1,y1]]*hwlist[{a1, a2, a3, a4, a6},P[x2,y2]]==0., 0., ArcCos[ heightpairing[{a1, a2, a3, a4, a6},P[x1,y1],P[x2,y2]]/ (Sqrt[hwlist[{a1, a2, a3, a4, a6},P[x1,y1]]]*Sqrt[hwlist[{a1, a2, a3, a4, a6},P[x2,y2]]]) ] ]; HeightRegulator[listofpoints_List]:= Det[Outer[HeightPairing,listofpoints,listofpoints]]/;Length[listofpoints]>0; HeightRegulator[listofpoints_List]:=1; HeightRegulator[]:=1; HeightRegulator[somepoints__]:= Det[ Outer[ HeightPairing,{somepoints},{somepoints}] ]; heightregulator[{a1_,a2_,a3_,a4,_a6_},listofpoints_List]:= Det[ Outer[ heightpairing[{a1, a2, a3, a4, a6},#1,#2]&, listofpoints,listofpoints ] ]/;Length[listofpoints]>0; heightregulator[{a1_,a2_,a3_,a4,_a6_},listofpoints_List]:=1; heightregulator[{a1_,a2_,a3_,a4,_a6_},{}]:=1; heightregulator[{a1_,a2_,a3_,a4,_a6_}]:=1; heightregulator[{a1_,a2_,a3_,a4,_a6_},somepoints__P]:= Det[ Outer[ heightpairing[{a1, a2, a3, a4, a6},#1,#2]&, {somepoints}, {somepoints} ] ]; HeightMatrix[]:={}; HeightMatrix[listofpoints__]:= MatrixForm[Outer[ HeightPairing,{listofpoints},{listofpoints}]]; HeightMatrix[listofpoints_List]:= MatrixForm[Outer[ HeightPairing,listofpoints,listofpoints]]; heightmatrix[{a1_,a2_,a3_,a4,_a6_}]:={}; heightmatrix[{a1_,a2_,a3_,a4,_a6_},somepoints__P]:= MatrixForm[ Outer[ heightpairing[{a1, a2, a3, a4, a6},#1,#2]&,{somepoints},{somepoints} ] ]; heightmatrix[{a1_,a2_,a3_,a4,_a6_},listofpoints_List]:= MatrixForm[ Outer[ heightpairing[{a1, a2, a3, a4, a6},#1,#2]&,listofpoints,listofpoints ] ]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Utilities for computing the reduction type :[font = input; initialization; preserveAspect; endGroup; ] *) changeacoeffs[{a1_,a2_,a3_,a4_,a6_},{r_,s_,t_}]:= {a1 + 2s, a2 - s a1 + 3r - s^2, a3 + r*a1 + 2t, a4 - s*a3 + 2r*a2 - (t+r*s)*a1 + 3r^2 - 2s*t, a6 + r*a4 + r^2*a2 + r^3 -t*a3 -t^2 - r*t*a1} numberofrootsink3[p_,d2_,d1_,d0_]:= Count[Mod[Map[(#^3+d2*#^2+d1*#+d0)&,Range[0,p-1]],p],0]/;PrimeQ[p]; numberofrootsink2[p_,d2_:1,d1_,d0_]:= Module[{ct=0}, If[p==2, If[Mod[d0,2]==0, ct++]; If[Mod[d0+d1+d2,2]==0, ct++]; ct, 1+JacobiSymbol[d1*d1-4d0*d2,p] ] ]/;PrimeQ[p]; numberofrootsink3[p_,d2_,d1_,d0_]:= Count[Mod[Map[(#^3+d2*#^2+d1*#+d0)&,Range[0,p-1]],p],0]/;PrimeQ[p]; discriminantcubic[a_,b_,c_]:= (-4a^3*c + a^2*b^2 + 18a*b*c -4b^3-27c^2); kodairatype[i_,sb_,sp_] := StringForm["``````",i,Subscript[sb], Superscript[sp]]; SINGULAR="Singular"; GOOD="Good"; SPLIT="Split"; NONSPLIT="Non-Split"; ADDITIVE="Additive"; NONMINIMAL="Non-minimal"; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Computing the Reduction Types :[font = input; initialization; preserveAspect; endGroup; ] *) Reductiontype[p_]:= ReductiontypeList[p][[7]]; reductiontype[p_]:= Reductiontype[p]; reductiontype[{a1_,a2_,a3_,a4_,a6_},p_]:= reductiontypewlist[{a1,a2,a3,a4,a6},p][[7]]; ReductiontypeList[p_]:= reductiontypewlist[{a[1],a[2],a[3],a[4],a[6]},p] reductiontypelist[p_]:= ReductiontypeList[p]; reductiontypelist[{a1_,a2_,a3_,a4_,a6_},p_]:= reductiontypewlist[{a1,a2,a3,a4,a6},p]; Clear[reductiontypewlist]; reductiontypewlist[{aa1_,aa2_,aa3_,aa4_,aa6_},p_]:= reductiontypewlist[{aa1,aa2,aa3,aa4,aa6},p]= Module[{a1,a2,a3,a4,a6,b2,b4,b6,b8,c4,c6,orddisc, m,n,rst,r0,s0,t0,f,c,disc}, {a1,a2,a3,a4,a6}={aa1,aa2,aa3,aa4,aa6}; b2:=a1^2+4a2; b4:=2a4+a1 a3; b6:=a3^2+4a6; b8:=a1^2 a6 + 4a2 a6 - a1 a3 a4 + a2 a3^2 - a4^2; c4:=b2^2-24b4; c6:=-b2^3 + 36 b2 b4 - 216 b6; disc:=(c4^3-c6^2)/1728; orddisc:= ord[disc,p]; aa[1,m_]:=a1/p^m; aa[2,m_]:=a2/p^m; aa[3,m_]:=a3/p^m; aa[4,m_]:=a4/p^m; aa[6,m_]:=a6/p^m; If[disc==0, Return[{p,orddisc,-1,"Sing",0,SINGULAR,-1}]]; (*Step1*) If[orddisc==0, Return[{p,0,0,kodairatype["I","0",""],1,GOOD,0}]]; (*Step2*) If[p==2, rst=Do[ If[ Mod[{a3+a1*r+2t, a4- s a3 +2r a2-(t+r s)a1 + 3r^2 - 2s t, a6+r a4 +r^2a2 + r^3 -t a3 -t^2 - r t a1},p ]=={0,0,0}, Return[{r,s,t}] ] ,{s,0,p-1},{r,0,p-1},{t,0,p-1}], rst=Do[ t=Mod[(-a3-r*a1)*(p+1)/2,p]; If[ Mod[{a3+a1*r+2t, a4- s a3 +2r a2-(t+r s)a1 + 3r^2 - 2s t, a6+r a4 +r^2a2 + r^3 -t a3 -t^2 - r t a1},p ]=={0,0,0}, Return[{r,s,t}] ] ,{s,0,p-1},{r,0,p-1}] ]; {a1,a2,a3,a4,a6}=changeacoeffs[{a1,a2,a3,a4,a6},rst]; If[Mod[b2,p]!=0, f=1; If[numberofrootsink2[p,a1,-a2]>0, c=orddisc;rst=SPLIT, c=Mod[orddisc-1,2]+1;rst=NONSPLIT]; Return[{p,orddisc,f,kodairatype["I",orddisc,""],c,rst, orddisc+10}] ]; (*Step 3*) If[Mod[a6,p^2]!=0, f=orddisc; c=1; Return[{p,orddisc,f,"II",c,ADDITIVE,2}]]; (*Step 4*) If[Mod[b8,p^3]!=0, f=orddisc-1; c=2; Return[{p,orddisc,f,"III",c,ADDITIVE,3}]]; (*Step 5*) If[Mod[b6,p^3]!=0, f=orddisc-2; c=If[numberofrootsink2[p,aa[3,1],-aa[6,2]]>0,3,1]; (*1+numberofrootsink2[p,aa[3,1],-aa[6,2]]; *) Return[{p,orddisc,f,"IV",c,ADDITIVE,4}]]; (*Step 6*) If[p!=3, Which[ p==2, r0=0;s0=Mod[a2,2];t0=Mod[a3/2,2], p==3, s0=Mod[a1,3]; (* Doesn't happen! *), True, s0=Mod[-a1*(p+1)/2,p]; r0=mod[(s0^2+s0*a1-a2)/3,p];t0=mod[(-a3-r0*a1)/2,p^2] ]; Do[If[ Mod[spare=changeacoeffs[{a1,a2,a3,a4,a6},{r,s,t}], {p,p,p^2,p^2,p^3} ]=={0,0,0,0,0}, {a1,a2,a3,a4,a6}=spare;Break[] ] ,{r,r0,p^2-1,p},{s,s0,p^2-1,p},{t,t0,p^2-1,p}], Do[If[ Mod[spare=changeacoeffs[{a1,a2,a3,a4,a6},{r,s,t}], {p,p,p^2,p^2,p^3} ]=={0,0,0,0,0}, {a1,a2,a3,a4,a6}=spare;Break[] ] ,{r,0,p^2-1},{s,0,p^2-1},{t,0,p^2-1}]]; If[Mod[discriminantcubic[aa[2,1],aa[4,2], aa[6,3]],p]!=0, f=orddisc-4;c=1+numberofrootsink3[p,aa[2,1],aa[4,2], aa[6,3]]; Return[{p,orddisc,f,kodairatype["I",0,"*"],c,ADDITIVE,-10}]; ]; (*Step 7*) If[((p==3) &&(Mod[aa[2,1],3]!=0 || Mod[aa[4,2],3]!=0)) || ((p!=3) && ((mod[aa[2,1]^2/3,p]!= Mod[aa[4,2],p]) || (mod[aa[2,1]^3/27,p]!= Mod[aa[6,3],p]))), (* We have a double root and a single root *) If[p!=3, r=p*Do[ If[Mod[3t^2+2aa[2,1]*t+aa[4,2],p]==0, If[Mod[t^3+aa[2,1]*t^2+aa[4,2]*t+aa[6,3],p]==0, Return[t], Return[2aa[2,1]-t] ] ] ,{t,0,p-1}], Which[ True, r=mod[aa[4,2]/aa[2,1],3]*3, (* If it gets this far, a[2,1]!=0 mod 3 *) Mod[aa[2,1],3]==Mod[aa[4,2],3]==0, r=0, Mod[aa[2,1],3]==Mod[aa[4,2],3]==Mod[-aa[6,3]-1,3], r=2p, True, r=p ] ]; {a1,a2,a3,a4,a6}=changeacoeffs[{a1,a2,a3,a4,a6},{r,0,0}]; m=2;n=2;nu=1;ppower=p*p; While[True, If[Mod[aa[3,m]^2+4*aa[6,2n],p]!=0, c=2+numberofrootsink2[p,aa[3,m],-aa[6,2n]];Break[] ]; If[p==2, t=aa[6,2n]*ppower, t=Mod[-aa[3,m]*(p+1)/2,p]*ppower]; {a1,a2,a3,a4,a6}=changeacoeffs[{a1,a2,a3,a4,a6},{0,0,t}]; nu++;m++; If[Mod[aa[4,n+1]^2-4*aa[2,1]*aa[6,2n+1],p]!=0, c=2+numberofrootsink2[p,mod[aa[4,n+1]/aa[2,1],p],mod[aa[6,2n+1]/aa[2,1],p]]; Break[] ]; If[p==2,r=aa[6,2n+1]*ppower, r=Mod[-aa[4,n+1]*(p+1)/2,p]*ppower]; {a1,a2,a3,a4,a6}=changeacoeffs[{a1,a2,a3,a4,a6},{r,0,0}]; nu++;n++;ppower*=p; ]; f=orddisc-4-nu; Return[{p,orddisc,f,kodairatype["I",nu,"*"],c,ADDITIVE,-nu-10}] ]; (*Step 8*) rst=Evaluate[{If[p==3,-a6/9,-a2*PowerMod[3,-1,p]],0,0}]; {a1,a2,a3,a4,a6}=changeacoeffs[{a1,a2,a3,a4,a6},rst]; If[Mod[aa[3,2]^2+4*aa[6,4],p]!=0, f=orddisc-6;c=1+numberofrootsink2[p,aa[3,2], -aa[6,4]]; Return[{p,orddisc,f,kodairatype["IV","","*"],c,ADDITIVE,-4}] ]; (*Step 9*) rst=Evaluate[{0,0,If[p==2,Mod[aa[6,4],2]*4,-aa[3,2]*(p+1)/2*p^2]}]; {a1,a2,a3,a4,a6}=changeacoeffs[{a1,a2,a3,a4,a6},rst]; If[Mod[aa[4,3],p]!=0, f=orddisc-7;c=2; Return[{p,orddisc,f,kodairatype["III","","*"],c,ADDITIVE,-3}] ]; (*Step 10*) If[ord[a6,p]<6, f=orddisc-8;c=1; Return[{p,orddisc,f,kodairatype["II","","*"],c,ADDITIVE,-2}] ]; (*Step 11*) (* Shouldn't happen if equation is minimal *) rst=Evaluate[{aa[1,1],aa[2,2],aa[3,3],aa[4,4],aa[6,6]}]; Message[ReductiontypeList::nonminimal,p,displayellipticcurve[rst]]; (*reductiontypewlist[p,rst]*) Return[{p, orddisc, f=1, "Not Min.",c=0, NONMINIMAL, 1}]; ] (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Computing the Conductor :[font = input; initialization; dontPreserveAspect; endGroup; ] *) ConductorList[] := TableForm[ Join[{{ " prime ","ord disc","ord cond","Kodaira","red type"}}, Map[ReductiontypeList[#][[{1,2,3,4,6}]]&, Select[Thread[DiscriminantFactored][[1]],Positive] ] ], TableAlignments->Center,TableSpacing->{0,2} ]; ConductorListFull[] := TableForm[ Join[{{ " prime ","ord disc","ord cond","c_p","Kodaira","red type"}}, Map[ReductiontypeList[#][[{1,2,3,5,4,6}]]&, Select[Thread[DiscriminantFactored][[1]],Positive] ] ], TableAlignments->Center,TableSpacing->{0,2} ] ; Conductor[] := conductorwlist[{a[1],a[2],a[3],a[4],a[6]}]; conductorlist[{a1_,a2_,a3_,a4_,a6_}] := TableForm[ Join[{{ " prime ","ord disc","ord cond","Kodaira","red type"}}, Map[reductiontypelist[{a1,a2,a3,a4,a6},#][[{1,2,3,4,6}]]&, Select[Thread[discriminantfactored[{a1,a2,a3,a4,a6}]][[1]],Positive] ] ], TableAlignments->Center,TableSpacing->{0,2} ] ; conductorlistfull[{a1_,a2_,a3_,a4_,a6_}] := TableForm[ Join[{{ " prime ","ord disc","ord cond","c_p","Kodaira","red type"}}, Map[reductiontypelist[{a1,a2,a3,a4,a6},#][[{1,2,3,5,4,6}]]&, Select[Thread[discriminantfactored[{a1,a2,a3,a4,a6}]][[1]],Positive] ] ], TableAlignments->Center,TableSpacing->{0,2} ] ; conductor[{a1_,a2_,a3_,a4_,a6_}] := conductorwlist[{a1,a2,a3,a4,a6}]; conductorwlist[{a1_, a2_, a3_, a4_, a6_}] := conductorwlist[{a1, a2, a3, a4, a6}] = (* Store Value *) Apply[Times, Map[ {#,ReductiontypeList[#][[3]]}&, Select[ Thread[discriminantfactored[{a1, a2, a3, a4, a6}]][[1]], Positive] ] /. {p_Integer,e_Integer} -> p^e ] SingularPoint[p_Integer] := Which[ Not[IntegerQ[Discriminant/p]], Origin, p==2, Switch[Mod[a[1],2], 0,P[Mod[a[4],2],Mod[a[2]*a[4]+a[6],2]], 1,P[Mod[a[3],2],Mod[a[3]+a[4],2]] ], p==3, Switch[Mod[b[2],3], 0,Mod[-b[6],3], 1,Mod[-b[4],3], 2,Mod[b[4],3] ] // P[#,Mod[a[1] #+a[3],3]]&, True, (* p >= 5 *) If[IntegerQ[c[4]/p], -PowerMod[12,-1,p]*b[2], -PowerMod[c[4],-1,p]*(b[2]*b[4]-18b[6]) ] //P[Mod[#,p],Mod[-(#a[1]+a[3])*(p+1)/2,p]]& ] singularpoint[{a1_,a2_,a3_,a4_,a6_},p_Integer] := Module[{b2,b4,b6,b8,c4,c6}, Which[ Not[IntegerQ[discfromaalist[{a1,a2,a3,a4,a6}]/p]], Origin, p==2, Switch[Mod[a1,2], 0,P[Mod[a4,2],Mod[a2*a4+a6,2]], 1,P[Mod[a3,2],Mod[a3+a4,2]] ], p==3, {b2,b4,b6,b8}=bbfromaalist[{a1,a2,a3,a4,a6}]; Switch[Mod[b2,3], 0,Mod[-b6,3], 1,Mod[-b4,3], 2,Mod[b4,3] ] // P[#,Mod[a1 #+a3,3]]&, True, (* p >= 5 *) {b2,b4,b6,b8}=bbfromaalist[{a1,a2,a3,a4,a6}]; {c4,c6}=ccfromaalist[{b2,b4,b6,b8}]; If[IntegerQ[c4/p], -PowerMod[12,-1,p]*b2, -PowerMod[c4,-1,p]*(b2*b4-18b6) ] //P[Mod[#,p],Mod[-(#a1+a3)*(p+1)/2,p]]& ] ] ChangeCoordinates[u_,r_,s_,t_] := ( (* DOES NOT CHANGE b's, c's, Discriminant *) Unprotect[a]; a[6] = (a[6]+r*(a[4]+r*(a[2]+r)) -t(a[3]+t+r*a[1]))/u^6; a[4] = (a[4]-s*a[3]+2r*a[2] -(t+r*s)*a[1]+3r^2-2s*t)/u^4; a[3] = (a[3]+r*a[1]+2t)/u^3; a[2] = (a[2]-s*a[1]+3r-s^2)/u^2; a[1] = (a[1]+2s)/u; Protect[a]; ); changecoordinates[{a1,a2,a3,a4,a6},u_,r_,s_,t_] := {(a1+2s)/u, (a2-s*a1+3r-s^2)/u^2, (a3+r*a1+2t)/u^3, (a4-s*a3+2r*a2 -(t+r*s)*a1+3r^2-2s*t)/u^4, (a6+r*(a4+r*(a2+r))-t(a3+t+r*a1))/u^6 } ChangeCoordinatesFull[u_,r_,s_,t_] := Module[{a1,a2,a3,a4,a6,myrules}, (* DOES CHANGE b's, c's, Discriminant *) Unprotect[a,b,c,Discriminant,DiscriminantFactored]; a[6] = (a[6]+r*(a[4]+r*(a[2]+r)) -t(a[3]+t+r*a[1]))/u^6; a[4] = (a[4]-s*a[3]+2r*a[2] -(t+r*s)*a[1]+3r^2-2s*t)/u^4; a[3] = (a[3]+r*a[1]+2t)/u^3; a[2] = (a[2]-s*a[1]+3r-s^2)/u^2; a[1] = (a[1]+2s)/u; a1=a[1]; a2=a[2]; a3=a[3]; a4=a[4]; a6=a[6]; b[2]=a1^2+4*a2; b[4]=a1*a3+2*a4; b[6]=a3^2+4*a6; b[8]=a1^2*a6+4*a2*a6-a1*a3*a4+a2*a3^2-a4^2; c[4]=c[4]/u^4(*b[2]^2-24*b[4]*); c[6]=c[6]/u^6(*-b[2]^3+36*b[2]*b[4]-216*b[6]*); Discriminant=Discriminant/u^12(*(c[4]^3-c[6]^2)/1728*); If[DiscriminantFactored!=Null, myrules=Map[{#[[1]],e_}->{#[[1]],e-12#[[2]]}&,FactorInteger[u]]; DiscriminantFactored=Select[DiscriminantFactored/.myrules,(#[[2]]!=0)&]; , ];(* Saves factoring Discriminant again *) Protect[a,b,c,Discriminant,DiscriminantFactored]; ]/;If[u==0,False,True,True] (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] NumberTheory Functions :[font = input; initialization; preserveAspect; endGroup; ] *) KroneckerSymbol[d_Integer?FundamentalDiscriminantQ,1]:=1; KroneckerSymbol[d_Integer,n_Integer]:= kroneckersymbol[d,n]/;FundamentalDiscriminantQ[d]&&(GCD[n,d]==1) KroneckerSymbol[d_Integer,n_Integer]:=0/;FundamentalDiscriminantQ[d]&&(GCD[n,d]!=1) KroneckerSymbol[d_Integer,n_Integer]:=Message[KroneckerSymbol::notcomp,d]; kroneckersymbol[d_Integer,1]:=1; kroneckersymbol[d_Integer,n_Integer]:= Block[{factors}, (*If[n<0, n=-Abs[d]*Ceiling[Abs[n/d]]+n, ];*) nn=Mod[n,d]+If[d<0,Abs[d],0]; If[GCD[d,n]>1, 0, factors=FactorInteger[n]; Which[ n==1, 1, Length[factors]>1, Apply[Times,factors/.{{p_Integer,e_Integer}->If[OddQ[e],kroneckersymbol[d,p],1]}], factors[[1,1]]==2, If[OddQ[factors[[1,2]]], JacobiSymbol[2,Abs[d]],1], True, If[OddQ[factors[[1,2]]], JacobiSymbol[d,factors[[1,1]]],1] ] ] ]; FundamentalDiscriminantQ[d_]:= Or[SquareFreeQ[d]&&(Mod[d,4]==1), (Mod[d,4]==0) && SquareFreeQ[d/4] && MemberQ[{2,3},Mod[d/4,4]] ]; FundamentalDiscriminant[d_Integer]:= Module[{z}, z=Apply[Times,FactorInteger[d]/.{a_Integer,b_Integer}->a^Mod[b,2]]; If[MemberQ[{2,3}, Mod[z,4]], 4z, z] ]; If[!ValueQ[SquareFreeQ], SquareFreeQ[d_Integer]:=(Max[FactorInteger[d]/.{{a_Integer,b_Integer}->b}]<2) , ]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Quadratic Twists :[font = input; initialization; preserveAspect; endGroup; ] *) QuadraticTwist[d_Integer]:= EllipticCurve[quadratictwist[astack[curvenumber],d]]/;FundamentalDiscriminantQ[d]; QuadraticTwist[d_Integer,n_Integer:Automatic]:= If[IntegerQ[n] && (1 <= n <= curvenumber), EllipticCurve[quadratictwist[d,astack[n]]], EllipticCurve[quadratictwist[d,astack[curvenumber]]] ]/;FundamentalDiscriminantQ[d]; QuadraticTwist[{a1_,a2_,a3_,a4_,a6_},d_Integer]:= EllipticCurve[quadratictwist[{a1,a2,a3,a4,a6},d]]/;FundamentalDiscriminantQ[d]; QuadraticTwist[d_Integer,___]:=Message[QuadraticTwist::notfunddiscrimant,d]; quadratictwist[{a1_,a2_,a3_,a4_,a6_},d_]:= Module[{nn}, If[a1==a3==0, nn={0,d*a2,0,d^2*a4,d^3*a6}, nn=ccfromaalist[{a1,a2,a3,a4,a6}]; nn={0,0,0,-27*nn[[1]]*d^2, -54nn[[2]]*d^3} ]; First[minimalequationandurst[nn]] ]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Modp :[font = input; initialization; dontPreserveAspect; endGroup; ] *) FrobeniusTrace[p_Integer] := frobeniustracewlist[{a[1],a[2],a[3],a[4],a[6]},p]; frobeniustrace[{a1_,a2_,a3_,a4_,a6_},p_Integer] := frobeniustracewlist[{a1,a2,a3,a4,a6},p]; frobeniustracewlist[{a1_,a2_,a3_,a4_,a6_},p_] := frobeniustracewlist[{a1,a2,a3,a4,a6},p] = (* Store Value *) Module[ {x,bb}, bb=bbfromaalist[{a1,a2,a3,a4,a6}]; Which[ p==2, 2 - Switch[{Mod[a3,2],Mod[a6,2]}, {1,1},0, {1,0},2, {0,_},1 ] - Switch[ {Mod[a1+a3,2],Mod[a2+a4+a6,2]}, {1,0},0, {1,1},2, {0,_},1 ], True, -Sum[JacobiSymbol[4x^3+bb[[1]]x^2+2bb[[2]]x+bb[[3]],p], {x,0,p-1} ] ] ]; SquareRootArray[p_]:= SquareRootArray[p]= Module[{indx,arr}, arr=Table[-1,{p}]; Do[arr[[Mod[indx*indx,p]]]=indx, {indx,1,Floor[p/2]}]; arr ](*/; PrimeQ[p]*); PointsModp[p_Integer]:= If[p==2, pointsmodp[{a[1],a[2],a[3],a[4],a[6]},2], If[PrimeQ[p], pointsmodp[{a[1],a[2],a[3],a[4],a[6]},p], "p must be prime"] ]; pointsmodp[{a1_,a2_,a3_,a4_,a6_},2]:= Module[{lst={Origin},z,k,kk}, Do[ If[Mod[k*k+a[1]*k*kk+a3*k-((kk+a2)*kk+a4)*kk+a6,2]==0, lst=Union[lst,{{kk,k}}] ,] ,{kk,0,1},{k,0,1}]; lst ]; pointsmodp[{a1_,a2_,a3_,a4_,a6_},p_Integer]:= Module[{lst={Origin},sqroot,z,k,temp,ff}, sqroot=SquareRootArray[p]; ff[xx_]:=((xx+a2)*xx+a4)*xx+a6; Do[ temp=a1*k+a3; z=Mod[temp*temp+4*ff[k],p]; If[z==0, AppendTo[lst,{k,Mod[-temp*(p+1)/2,p]}], If[(z=sqroot[[z]])==-1, , lst=Join[lst, {{k,Mod[(-temp+z)*(p+1)/2,p]}, {k,Mod[(-temp-z)*(p+1)/2,p]}} ] ] ] ,{k,0,p-1} ]; lst ]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Periods :[font = input; initialization; dontPreserveAspect; endGroup; ] *) RealPeriod[] := realperiodwlist[{a[1],a[2],a[3],a[4],a[6]}]; realperiodwlist[{a1_, a2_, a3_, a4_, a6_}] := realperiodwlist[{a1, a2, a3, a4, a6}] = (* Save value *) Module[{rootlist,realroots,e1,e2,e3,c4,c6}, {c4,c6}=ccfromaalist[{a1, a2, a3, a4, a6}]; rootlist = (x/.{ToRules[ NRoots[x^3-27c4*x-54c6==0,x]]})/36; realroots = Sort[Select[rootlist,Im[#]==0&]]; Which[ Length[realroots]==3, ( {e1,e2,e3} = realroots; 2N[Pi/ArithmeticGeometricMean[ Sqrt[e3-e1],Sqrt[e3-e2]]] ), Length[realroots]==1, ( e1 = realroots[[1]]; 2N[Pi/ArithmeticGeometricMean[ 2(3e1^2-c4/48)^(.25), Sqrt[3e1+2Sqrt[3e1^2-c4/48]]]]) ] ]; realperiod[{a1_, a2_, a3_, a4_, a6_}]:=realperiodwlist[{a1,a2,a3,a4,a6}] ImaginaryPeriod[] := imaginaryperiodwlist[{a[1],a[2],a[3],a[4],a[6]}]/;(Discriminant>0); imaginaryperiod[{a1_, a2_, a3_, a4_, a6_}]:= imaginaryperiodwlist[{a1,a2,a3,a4,a6}]/;(discriminant[{a1,a2,a3,a4,a6}]>0) ImaginaryPeriod::noncomp = "The Discriminant must be positive to compute the ImaginaryPeriod."; ImaginaryPeriod[]:=Message[ImaginaryPeriod::noncomp]; imaginaryperiodwlist[{a1_, a2_, a3_, a4_, a6_}] := imaginaryperiodwlist[{a1, a2, a3, a4, a6}] = (* Save value *) Module[{rootlist,realroots,e1,e2,e3,c4,c6}, {c4,c6}=ccfromaalist[{a1, a2, a3, a4, a6}]; rootlist = (x/.{ToRules[ NRoots[x^3-27c4*x-54c6==0,x]]})/36; realroots = Sort[Select[rootlist,Im[#]==0&]]; Which[ Length[realroots]==3, ( {e1,e2,e3} = realroots; I*N[Pi,20]/ArithmeticGeometricMean[ Sqrt[e3-e1],Sqrt[e2-e1]] ), Length[realroots]==1, ( e1 = realroots[[1]]; I*N[Pi,20]/ArithmeticGeometricMean[ 2(3e1^2-c4/48)^(.25), Sqrt[3e1+2Sqrt[3e1^2-c4/48]]]) ] ]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] L-Series :[font = input; initialization; dontPreserveAspect; endGroup; ] *) FudgeFactor[]:= Apply[ Times, Thread[ Map[ReductiontypeList, Select[ Thread[DiscriminantFactored][[1]], #>1&] ] ][[5]] ]; fudgefactor[{a1_, a2_, a3_, a4_, a6_}]:= Apply[ Times, Thread[ Map[reductiontypewlist[{a1, a2, a3, a4, a6} ,#]&, Select[ Thread[discriminantfactored[{a1, a2, a3, a4, a6}]][[1]], #>1&] ] ][[5]] ]; LSeriesCoefficient[1] = 1; LSeriesCoefficient[n_Integer] := lseriescoefficientwlist[{a[1],a[2],a[3],a[4],a[6]},n]; lseriescoefficient[{a1_,a2_,a3_,a4_,a6_},n_] := lseriescoefficientwlist[{a1,a2,a3,a4,a6},n]; lseriescoefficientwlist[{a1_,a2_,a3_,a4_,a6_},1] :=1; lseriescoefficientwlist[{a1_,a2_,a3_,a4_,a6_},n_] := lseriescoefficientwlist[{a1,a2,a3,a4,a6},n] = (* Store Value *) Module[{factors}, factors = FactorInteger[n]; Which[ Length[factors] > 1, (* n not a prime power *) Apply[Times, Map[lseriescoefficientwlist[{a1,a2,a3,a4,a6},#]&, factors /. {p_Integer,e_Integer}->p^e] ], factors[[1,2]] > 1, (* n is a prime power *) If[IntegerQ[discriminant[{a1,a2,a3,a4,a6}]/factors[[1,1]] ], frobeniustracewlist[{a1,a2,a3,a4,a6},factors[[1,1]]]^factors[[1,2]], (lseriescoefficientwlist[{a1,a2,a3,a4,a6},#1^(#2-1)]* lseriescoefficientwlist[{a1,a2,a3,a4,a6},#1] - #1* lseriescoefficientwlist[{a1,a2,a3,a4,a6},#1^(#2-2)])& [factors[[1,1]],factors[[1,2]] ] ], True, (* n is a prime *) frobeniustracewlist[{a1,a2,a3,a4,a6},factors[[1,1]]] ] ]; LSeries[1,opts___] := Block[ {NumberOfTerms}, Module[{n,ratio,conductor=Conductor[]}, NumberOfTerms = NumberOfTerms /. {opts} /. Options[LSeries]; If[FunctionalEquationSign[NumberOfTerms] == -1, Return[0] ]; ratio = Exp[-2*Pi/Sqrt[conductor]] //N ; 2*Sum[ LSeriesCoefficient[n]*ratio^n/n,{n,1,NumberOfTerms} ] ] ]; lseries[{a1_,a2_,a3_,a4_,a6_},1,opts___] := Block[ {NumberOfTerms}, Module[{n,ratio,conductor=conductorwlist[{a1,a2,a3,a4,a6}]}, NumberOfTerms = NumberOfTerms /. {opts} /. Options[LSeries]; If[functionalequationsign[{a1,a2,a3,a4,a6},NumberOfTerms] == -1, Return[0] ]; ratio = Exp[-2*Pi/Sqrt[conductor]] //N ; 2*Sum[ lseriescoefficientwlist[{a1,a2,a3,a4,a6},n]*ratio^n/n,{n,1,NumberOfTerms} ] ] ]; FunctionalEquationSign[NumberOfTerms_:20]:= If[functionalequationsign$0=!=0, functionalequationsign$0, functionalequationsign$0=functionalequationsign[{a[1],a[2],a[3],a[4],a[6]},NumberOfTerms]; If[functionalequationsign$0=0, Message[FunctionalEquationSign::notcomp,NumberOfTerms]]; functionalequationsign$0 ]; FunctionalEquationSign[opts__Rule|opts__RuleDelayed]:= Block[ {NumberOfTerms}, NumberOfTerms = NumberOfTerms/.{opts}; If[NumberQ[NumberOfTerms], Return[FunctionalEquationSign[NumberOfTerms]], Return[FunctionalEquationSign[]] ] ]; functionalequationsign[{a1_,a2_,a3_,a4_,a6_},NumberOfTerms_:20]:= Module[{singulartraces,alpha,q,n,sum,w}, (* Check if only multiplicative reduction *) singulartraces = Map[reductiontypewlist[{a1,a2,a3,a4,a6},#][[6]]&, Select[ Thread[discriminantfactored[{a1,a2,a3,a4,a6}] ][[1]],#>1&] ]; If[ (w=Count[singulartraces,SPLIT])+ (Count[singulartraces,NONSPLIT]) ==Length[singulartraces], Return[(-1)^(1+w)] ]; (* Some additive reduction *) alpha = .5 + .5*Sqrt[3]*I; q = Exp[2*Pi*I*alpha/Sqrt[conductor[{a1,a2,a3,a4,a6}]]] //N; sum=Sum[lseriescoefficientwlist[{a1,a2,a3,a4,a6},n]*q^n, {n,1,NumberOfTerms}] //N; w = Conjugate[sum]/(sum*alpha^2)//N; (* Note sign for modular form is opposite to sign for functional equation *) Which[ Abs[Re[w] - 1] < 10^-3, Return[-1], Abs[Re[w] + 1] < 10^-3, Return[1], True, Return[0] ] ]; FunctionalEquationSignTwist[eps_Integer,Cond_Integer,d_Integer,opts___]:= kroneckersymbol[d,Cond]*eps*Sign[d]/;(GCD[d,Cond]==1)&&FundamentalDiscriminantQ[d] FunctionalEquationSignTwist[eps_Integer,Cond_Integer,d_Integer,opts___]:= Block[{epsilon,nn=curvenumber}, QuadraticTwist[d]; epsilon=FunctionalEquationSign[opts]; OldCurve[nn]; epsilon ]/;((GCD[d,Cond]!=1)&&FundamentalDiscriminantQ[d]) FunctionalEquationSignTwist[{a1_,a2_,a3_,a4_,a6_},d_Integer,opts___]:= Block[{epsilon,nn=curvenumber}, quadratictwist[{a1,a2,a3,a4,a6},d]; epsilon=FunctionalEquationSign[opts]; OldCurve[nn]; epsilon ]/;((GCD[d,Cond]!=1)&&FundamentalDiscriminantQ[d]) functionalequationsigntwist[{a1_,a2_,a3_,a4_,a6_},d_Integer,opts___]:= functionalequationsign[quadratictwist[{a1,a2,a3,a4,a6},d] ]/;((GCD[d,Cond]!=1)&&FundamentalDiscriminantQ[d]) FunctionalEquationSignTwist[a___]:= Message[FunctionalEquationSignTwist::notfundisc]; bsd[{a1_,a2_,a3_,a4_,a6_}]:= realperiodwlist[{a1,a2,a3,a4,a6}]*fudgefactor[{a1,a2,a3,a4,a6}]/Length[torsiongrouppoints[{a1,a2,a3,a4,a6}]]^2; bsd[{a1_,a2_,a3_,a4_,a6_},a__P]:= ({bsd[{a1,a2,a3,a4,a6}],heightregulator[{a1,a2,a3,a4,a6},a]*2^Length[{a}]}/. {{aa_,bb_}->aa*bb}) /;minimalQ[{a1,a2,a3,a4,a6}]; bsdnonminimal[{a1_,a2_,a3_,a4_,a6_}]:= bsd[First[minimalequationandurst[{a1,a2,a3,a4,a6}]]]/;!minimalQ[{a1,a2,a3,a4,a6}]; BSD::nonlist="BSD requires a list of points."; BSD::nonminimal="BSD requires that the curve be minimal"; BSD[]:=bsd[{a[1],a[2],a[3],a[4],a[6]}]/;MinimalQ[]; BSD[]:=Message[BSD::nonminimal]; BSD[a__P]:= ({bsd[{a[1],a[2],a[3],a[4],a[6]}],HeightRegulator[a]*2^Length[{a}]}/. {{aa_,bb_}->aa*bb}) /;MinimalQ[]; BSD[a_List]:= ({bsd[{a[1],a[2],a[3],a[4],a[6]}],HeightRegulator[a]*2^Length[a]}/. {{aa_,bb_}->aa*bb}) /; MinimalQ[]; BSD[a_List|a_P]:=Message[BSD::nonminimal]; BSD[__]:=Message[BSD::nonlist] (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Modular Forms :[font = input; initialization; dontPreserveAspect; endGroup; ] *) PullBackLSeries[nn_Integer]:= Apply[Plus,Map[(LSeriesCoefficient[#]*q^#)&,Range[nn]]]; PullBackLSeries[{a1,a2,a3,a4,a6},nn_Integer]:= Apply[Plus, Map[(lseriescoefficient[{a1,a2,a3,a4,a6},#]*q^#)&,Range[nn]] ]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Graphing :[font = input; initialization; dontPreserveAspect; endGroup; ] *) PlotEllipticCurve[xmax_:5,options___] := Module[{e1,e2,e3,x,SQRT,rootlist,realroots,leftpiece,rightpiece,label}, SQRT[s_] := Re[Sqrt[N[s]]]; label = PlotLabel /. {options} /. Options[PlotEllipticCurve]; rootlist = (x/.{ToRules[ NRoots[4x^3+b[2]x^2+2b[4]x+b[6]==0,x]]}); realroots = Sort[Select[rootlist,Im[#]==0&]]; If[realroots[[1]]==realroots[[2]], realroots={realroots[[3]]}]; Which[ Length[realroots]==3, ( {e1,e2,e3} = realroots; leftpiece = Plot[{ .5*(-a[1]x-a[3]+SQRT[4x^3+b[2]x^2+2b[4]x+b[6]]), .5*(-a[1]x-a[3]-SQRT[4x^3+b[2]x^2+2b[4]x+b[6]])}, {x,e1,e2},options,DisplayFunction->Identity]; rightpiece = Plot[{ .5*(-a[1]x-a[3]+SQRT[4x^3+b[2]x^2+2b[4]x+b[6]]), .5*(-a[1]x-a[3]-SQRT[4x^3+b[2]x^2+2b[4]x+b[6]])}, {x,e3,xmax},options,DisplayFunction->Identity]; Show[leftpiece,rightpiece, PlotLabel -> label, DisplayFunction -> $DisplayFunction]), Length[realroots]==1, ( e1 = realroots[[1]]; Plot[{ .5*(-a[1]x-a[3]+SQRT[4x^3+b[2]x^2+2b[4]x+b[6]]), .5*(-a[1]x-a[3]-SQRT[4x^3+b[2]x^2+2b[4]x+b[6]])}, {x,e1,xmax},options, PlotLabel -> label]) ] ]; plotellipticcurve[{a1_,a2_,a3_,a4_,a6_},xmax_:5,options___] := Module[{e1,e2,e3,x,SQRT,rootlist,realroots,leftpiece,rightpiece,label, b2,b4,b6,b8}, SQRT[s_] := Re[Sqrt[N[s]]]; label = PlotLabel /. {options} /. Options[PlotEllipticCurve]; {b2,b4,b6,b8}=bbfromaalist[{a1,a2,a3,a4,a6}]; rootlist = (x/.{ToRules[ NRoots[4x^3+b2x^2+2b4x+b6==0,x]]}); realroots = Sort[Select[rootlist,Im[#]==0&]]; If[realroots[[1]]==realroots[[2]], realroots={realroots[[3]]}]; Which[ Length[realroots]==3, ( {e1,e2,e3} = realroots; leftpiece = Plot[{ .5*(-a1x-a3+SQRT[4x^3+b2x^2+2b4x+b6]), .5*(-a1x-a3-SQRT[4x^3+b2x^2+2b4x+b6])}, {x,e1,e2},options,DisplayFunction->Identity]; rightpiece = Plot[{ .5*(-a1x-a3+SQRT[4x^3+b2x^2+2b4x+b6]), .5*(-a1x-a3-SQRT[4x^3+b2x^2+2b4x+b6])}, {x,e3,xmax},options,DisplayFunction->Identity]; Show[leftpiece,rightpiece, PlotLabel -> label, DisplayFunction -> $DisplayFunction]), Length[realroots]==1, ( e1 = realroots[[1]]; Plot[{ .5*(-a1x-a3+SQRT[4x^3+b2x^2+2b4x+b6]), .5*(-a1x-a3-SQRT[4x^3+b2x^2+2b4x+b6])}, {x,e1,xmax},options, PlotLabel -> label]) ] ]; GraphPoints[pointlist:{P[_,_]..},opts___] := Graphics[{ PointSize[0.013], Map[ {Point[{#[[1]],#[[2]]}], Switch[ PointLabels/.{opts}/.Options[GraphPoints], Automatic, Text[StringForm["P[``,``]",#[[1]],#[[2]]], {#[[1]],#[[2]]},{-1.4,0}], None, {}, _, Print["Bad point label option"] ] } &, pointlist ] } ]; graphpoints[a___]:=GraphPoints[a]; GraphTangents[pointlist:{P[_,_]..},opts___] := { Map[Block[{Slope,x,y,dx,dy}, x=#[[1]]; y=#[[2]]; dx=dx/.{opts}/.Options[GraphTangents]; dy=dy/.{opts}/.Options[GraphTangents]; If[2y+a[1]x+a[3]==0, Graphics[Line[{{x,y-dy},{x,y+dy}}]], (Slope=(3x^2+2a[2]x+a[4]-a[1]y)/(2y+a[1]x+a[3]); Graphics[ Line[{{x-dx,y-Slope*dx},{x+dx,y+Slope*dx}}]]) ] ]&,pointlist ] , If[DisplayPoints/.{opts}/.Options[GraphTangents], GraphPoints[pointlist,opts],{},{}] }; graphtangents[{a1_,a2_,a3_,a4_,a6_},pointlist:{P[_,_]..},opts___] := { Map[Block[{Slope,x,y,dx,dy}, x=#[[1]]; y=#[[2]]; dx=dx/.{opts}/.Options[GraphTangents]; dy=dy/.{opts}/.Options[GraphTangents]; If[2y+a1x+a3==0, Graphics[Line[{{x,y-dy},{x,y+dy}}]], (Slope=(3x^2+2a2x+a4-a1y)/(2y+a1x+a[3]); Graphics[ Line[{{x-dx,y-Slope*dx},{x+dx,y+Slope*dx}}]]) ] ]&,pointlist ] , If[DisplayPoints/.{opts}/.Options[GraphTangents], GraphPoints[pointlist,opts],{},{}] }; GraphSecants[pointlist:{{P[_,_],P[_,_]}..},opts___] := { Map[ Block[{x1,y1,x2,y2,dx,dy,dt}, x1=#[[1,1]]; y1=#[[1,2]]; x2=#[[2,1]]; y2=#[[2,2]]; dx=x2-x1; dy=y2-y1; dt=Extension/.{opts}/.Options[GraphSecants]; Graphics[Line[{{x1-dx dt,y1-dy dt},{x2+dx dt,y2+dy dt}}]] ]&, pointlist] , If[DisplayPoints/.{opts}/.Options[GraphTangents], GraphPoints[Flatten[pointlist],opts],{},{}] }; graphsecants[a___]:=GraphSecants[a]; (* :[font = text; inactive; initialization; Cclosed; preserveAspect; startGroup; ] Equations of Secants and Tangents and the curve iteself. :[font = input; initialization; preserveAspect; endGroup; ] *) TangentLine[P[xx_,yy_]]:= Module[{slop,intcpt,temp,line}, temp=(2yy+a[1]xx+a[3]); Switch[temp, 0, line=(x==xx), _, slop=(3xx^2+2a[2]xx+a[4]-a[1]yy)/temp; intcpt=yy-slop*xx; line=(y==slop*x+intcpt); ]; line ] tangentline[{a1_,a2_,a3_,a4_,a6_},P[xx_,yy_]]:= Module[{slop,intcpt,temp,line}, temp=(2yy+a1xx+a3); Switch[temp, 0, line=(x==xx), _, slop=(3xx^2+2a2*xx+a4-a1*yy)/temp; intcpt=yy-slop*xx; line=(y==slop*x+intcpt); ]; line ] SecantLine[P[xx_,yy_],P[zz_,ww_]]:= Module[{slop,intcpt,line}, If[zz==xx, If[yy==ww, line=TangentLine[P[xx,yy]], line=(x==xx)], slop=(ww-yy)/(zz-xx); intcpt=yy-slop*xx; line=(y==slop*x+intcpt); ]; line ] secantline[a___]:=SecantLine[a]; EquationCurve[xx_:x,yy_:y]:= (yy^2+a[1]xx*yy+a[3]*yy==xx^3+a[2]*xx^2+a[4]*xx+a[6]) (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Legendre Elliptic Curves y^2=x(x-1)(x-lambda) :[font = input; initialization; dontPreserveAspect; endGroup; ] *) OmegaOne[lambda_]:= Module[{lst={},lst2={},lambdaa}, If[(lambda==0)||(lambda==1), Message[OmegaOne::sing,lambda], lst={lambda,1/lambda,(lambda-1)/lambda, lambda/(lambda-1),1/(lambda-1),1-lambda}; lst2=Select[lst, ((Abs[#]<1)&&(Abs[#-1]<1))&]; If[Length[lst2]>=1, lambdaa=First[lst2]; I*Pi*Hypergeometric2F1[-0.5, -0.5, 1, 1-lambdaa]//N, Message[OmegaOne::notcomp,lambdaa] ] ] ]; OmegaTwo[lambda_]:= Module[{lst={},lst2={},lambdaa}, If[(lambda==0)||(lambda==1), Message[OmegaOne::sing,lambda], lst={lambda,1/lambda,(lambda-1)/lambda, lambda/(lambda-1),1/(lambda-1),1-lambda}; lst2=Select[lst, ((Abs[#]<1)&&(Abs[#-1]<1))&]; If[Length[lst2]>=1, lambdaa=First[lst2]; Pi*Hypergeometric2F1[-0.5, -0.5, 1, lambda]//N, Message[OmegaOne::notcomp,lambdaa] ] ] ]; EllipticCurveLegendre[lambda_,opts___]:=EllipticCurve[0,-lambda-1,0,lambda,0,opts]; (* :[font = text; inactive; initialization; Cclosed; dontPreserveAspect; startGroup; ] Display the a, b, and c coefficients, the Discriminant, and the j-invariant :[font = input; initialization; dontPreserveAspect; endGroup; endGroup; ] *) Displayacoeffs[]:= Block[{}, Print[" a coefficients "]; Map[Print["a",#," ",Factor[Simplify[a[#]]]]&,{1,2,3,4,6}]; ]; Displaybcoeffs[]:= Block[{}, Print[" b coefficients "]; Map[Print["b",#," ",Factor[Simplify[b[#]]]]&,{2,4,6,8}]; ]; Displayccoeffs[]:= Block[{}, Print[" c coefficients "]; Map[Print["c",#," ",Factor[Simplify[c[#]]]]&,{4,6}]; ]; DisplayDiscj[]:= Module[{zz}, Print[" Discriminant = ",Discriminant," = ", FactorInteger[Simplify[Expand[Discriminant]]] ]; Print[" j = ",zz=Simplify[Expand[c[4]^3/Discriminant]]," = ", FactorInteger[zz]]; ]; (* :[font = subsection; inactive; initialization; Cclosed; preserveAspect; startGroup; ] Epilogue :[font = input; initialization; preserveAspect; ] *) End[]; (* :[font = input; initialization; preserveAspect; endGroup; endGroup; ] *) EndPackage[]; (* ^*)