Ticket #12766: 12766.patch

File 12766.patch, 8.1 KB (added by roed, 9 years ago)
  • sage/schemes/elliptic_curves/ell_generic.py

    # HG changeset patch
    # User David Roe <roed.math@gmail.com>
    # Date 1332977584 -3600
    # Node ID 472b4ab00375b3804d3e9053739f937c6fc3e9bb
    # Parent  a5ec30a5f18a669fb792109fa8f69804928cedf4
    Improved plotting for elliptic curves
    
    diff --git a/sage/schemes/elliptic_curves/ell_generic.py b/sage/schemes/elliptic_curves/ell_generic.py
    a b  
    5151import sage.groups.additive_abelian.additive_abelian_group as groups
    5252import sage.groups.generic as generic
    5353import sage.plot.all as plot
     54from sage.plot.plot import generate_plot_points
    5455
    5556import sage.rings.arith as arith
    5657import sage.rings.all as rings
     
    24942495
    24952496   
    24962497    # Plotting
    2497    
    24982498
    2499     def plot(self, xmin=None, xmax=None, **args):
     2499
     2500    def plot(self, xmin=None, xmax=None, components='both', **args):
    25002501        """
    25012502        Draw a graph of this elliptic curve.
    25022503       
     
    25042505       
    25052506        -  ``xmin, xmax`` - (optional) points will be computed at
    25062507           least within this range, but possibly farther.
     2508
     2509        - ``components`` - a string, one of the following:
     2510
     2511          - ``both`` -- (default), scale so that both bounded and
     2512            unbounded components appear
     2513
     2514          - ``bounded`` -- scale the plot to show the bounded
     2515            component.  Raises an error if there is only one real
     2516            component.
    25072517       
     2518          - ``unbounded`` -- scale the plot to show the unbounded
     2519            component, including the two flex points.
     2520
    25082521        -  ``**args`` - all other options are passed to the
    25092522           line graphing primitive.
    25102523       
     
    25232536            RR._coerce_(K(1))
    25242537        except TypeError:
    25252538            raise NotImplementedError, "Plotting of curves over %s not implemented yet"%K
     2539        if components not in ['both', 'bounded', 'unbounded']:
     2540            raise ValueError("component must be one of 'both', 'bounded' or 'unbounded'")
     2541       
    25262542        a1, a2, a3, a4, a6 = self.ainvs()
    2527         R = rings.PolynomialRing(rings.RealField(), 'x')
    2528         x = R.gen()
    2529         d = 4*x**3 + (a1**2 + 4*a2)*x**2 + (2*a3*a1 + 4*a4)*x + (a3**2 + 4*a6)
     2543        d = self.division_polynomial(2)
    25302544        # Internal function for plotting first branch of the curve
    25312545        f1 = lambda z: (-(a1*z + a3) + sqrt(abs(d(z))))/2
    25322546        # Internal function for plotting second branch of the curve
    25332547        f2 = lambda z: (-(a1*z + a3) - sqrt(abs(d(z))))/2
    2534         r = d.roots(multiplicities=False)
     2548        r = d.roots(RR, multiplicities=False)
    25352549        r.sort()
    2536         if xmax is None:
    2537             xmax = r[-1] + 2
    2538         xmax = max(xmax, r[-1]+2)
    2539         if xmin is None:
    2540             xmin = r[0]  - 2
    2541         xmin = min(xmin, r[0]-2)
    2542         if len(r) == 1:
     2550        if components == 'bounded' and len(r) == 1:
     2551            raise ValueError("no bounded component for this curve")
     2552        if xmin is None or xmax is None:
     2553            xmins = []
     2554            xmaxs = []
     2555            if components in ['both','bounded'] and len(r) > 1:
     2556                xmins.append(r[0])
     2557                xmaxs.append(r[1])
     2558            # The following 3 is an aesthetic choice.
     2559            if components == 'unbounded' or components == 'both' and (len(r) == 1 or r[2] - r[1] > 3*(r[1] - r[0])):
     2560                flex = self.division_polynomial(3).roots(RR, multiplicities=False)
     2561                flex.sort()
     2562                flex = flex[-1]
     2563                xmins.append(r[-1])
     2564                # The doubling here is an aesthetic choice
     2565                xmaxs.append(flex + 2*(flex - r[-1]))
     2566            elif components == 'both':
     2567                # First the easy part.
     2568                xmins.append(r[-1])
     2569                # There are two components and the unbounded component
     2570                # is not too far from the bounded one.  We scale so
     2571                # that the unbounded component is twice as tall as the
     2572                # bounded component.  The y values corresponding to
     2573                # horizontal tangent lines are determined as follows.
     2574                # We implicitly differentiate the equation for this
     2575                # curve and get
     2576                # 2 yy' + a1 y + a1 xy' + a3 y' = 3 x^2 + 2a2 x + a4
     2577
     2578                R = RR['x']
     2579                x = R.gen()
     2580                if a1 == 0:
     2581                    # a horizontal tangent line can only occur at a root of
     2582                    Ederiv = 3*x**2 + 2*a2*x + a4
     2583                else:
     2584                    # y' = 0  ==>  y = (3*x^2 + 2*a2*x + a4) / a1
     2585                    y = (3*x**2 + 2*a2*x + a4) / a1
     2586                    Ederiv = y**2 + a1*x*y + a3*y - (x**3 + a2*x**2 + a4*x + a6)
     2587                critx = [a for a in Ederiv.roots(RR, multiplicities=False) if r[0] < a < r[1]]
     2588                if len(critx) == 0:
     2589                    raise RuntimeError("No horizontal tangent lines on bounded component")
     2590                # The 2.5 here is an aesthetic choice
     2591                ymax = 2.5 * max([f1(a) for a in critx])
     2592                ymin = 2.5 * min([f2(a) for a in critx])
     2593                top_branch = ymax**2 + a1*x*ymax + a3*ymax - (x**3 + a2*x**2 + a4*x + a6)
     2594                bottom_branch = ymin**2 + a1*x*ymin + a3*ymin - (x**3 + a2*x**2 + a4*x + a6)
     2595                xmaxs.append(max(top_branch.roots(RR,multiplicities=False) + bottom_branch.roots(RR,multiplicities=False)))
     2596            xmins = min(xmins)
     2597            xmaxs = max(xmaxs)
     2598            span = xmaxs - xmins
     2599            if xmin is None:
     2600                xmin = xmins - .02*span
     2601            if xmax is None:
     2602                xmax = xmaxs + .02*span
     2603        elif xmin >= xmax:
     2604            raise ValueError("xmin must be less than xmax")
     2605
     2606        I = []
     2607        if components in ['unbounded', 'both'] and xmax > r[-1]:
    25432608            # one real root; 1 component
    2544             I = [(r[0],xmax)]
    2545         else:
    2546             # three real roots; 2 components
    2547             I = [(r[0],r[1]), (r[2],xmax)]
    2548         I = [(max(a,xmin),min(b,xmax)) for a,b in I]
     2609            if xmin <= r[-1]:
     2610                I.append((r[-1],xmax,'<'))
     2611            else:
     2612                I.append((xmin, xmax,'='))
     2613        if components in ['bounded','both'] and len(r) > 1 and (xmin < r[1] or xmax > r[0]):
     2614            if xmin <= r[0]:
     2615                if xmax >= r[1]:
     2616                    I.append((r[0],r[1],'o'))
     2617                else:
     2618                    I.append((r[0],xmax,'<'))
     2619            elif xmax >= r[1]:
     2620                I.append((xmin, r[1], '>'))
     2621            else:
     2622                I.append((xmin, xmax, '='))
    25492623
    25502624        g = plot.Graphics()
    2551         try:
    2552             plot_points = int(args['plot_points'])
    2553             del args['plot_points']
    2554         except KeyError:
    2555             plot_points = 100
    2556            
     2625        plot_points = int(args.pop('plot_points',200))
     2626        adaptive_tolerance = args.pop('adaptive_tolerance',0.01)
     2627        adaptive_recursion = args.pop('adaptive_recursion',5)
     2628        randomize = args.pop('randomize',True)
    25572629        for j in range(len(I)):
    2558             a,b = I[j]
    2559             delta = (b-a)/float(plot_points)
    2560             v = []
    2561             w = []
    2562             for i in range(plot_points):
    2563                 x = a + delta*i
    2564                 v.append((x, f1(x)))
    2565                 w.append((x, f2(x)))
    2566             v.append((b,f1(b)))
    2567             w.append((b,f2(b)))
    2568             if len(I) == 2 and j == 0:  # two components -- the oh.
     2630            a,b,shape = I[j]
     2631            v = generate_plot_points(f1, (a, b), plot_points, adaptive_tolerance, adaptive_recursion, randomize)
     2632            w = generate_plot_points(f2, (a, b), plot_points, adaptive_tolerance, adaptive_recursion, randomize)
     2633            if shape == 'o':
    25692634                g += plot.line(v + list(reversed(w)) + [v[0]], **args)
     2635            elif shape == '<':
     2636                g += plot.line(list(reversed(v)) + w, **args)
     2637            elif shape == '>':
     2638                g += plot.line(v + list(reversed(w)), **args)
    25702639            else:
    2571                 g += plot.line(list(reversed(v)) + w, **args)
     2640                g += plot.line(v, **args)
     2641                g += plot.line(w, **args)
    25722642        return g
    25732643
    25742644    def formal_group(self):