# Algebraic curve plotter

An algebraic curve is a certain "1D subset" of the "2D-plane" that can be described as set of zeros {(x,y) in R^2 : f(x,y)=0 }of a polynomial f. Here we consider the 2D-plane as the real plane R^2 such that we can easily imagine what such a curve could look like, basically a thing you can draw with a pencil.

Examples:

• 0 = x^2 + y^2 -1 a circle of radius 1
• 0 = x^2 + 2y^2 -1 an ellipse
• 0 = xy a cross shape, basically the union of the x-axis and the y-axis
• 0 = y^2 - x a parabola
• 0 = y^2 - (x^3 - x + 1) an elliptic curve
• 0 = x^3 + y^3 - 3xy the folium of Descartes
• 0 = x^4 - (x^2 - y^2) a lemniscate
• 0 = (x^2 + y^2)^2 - (x^3 - 3xy^2) a trifolium
• 0 = (x^2 + y^2 - 1)^3 + 27x^2y^2 an astroid

Given a polynomial f (as defined below), and the x/y-ranges, output an black and white image of at least 100x100 pixels that shows the curve as black line on a white background.

### Details

Color: You can use any two other colours of your choice, it should just be easy to tell them apart.

Plot: Instead of a pixel image you can also output this image as ascii-art, where the background "pixels" should be space/underline or an other character that "looks empty" and the line can be made of a character that looks "full" like M or X or #.

You do not have to worry about aliasing.

You only need to plot lines where the sign of the polynomial changes from one side of the line to the other (that means you could e.g. use the marching square algorithm), you do not have to correctly plot "pathological cases like 0 = x^2 where the sign does not change when going from one side of the line to the other. But the line should be continuous and separating the regions of the different signs of f(x,y).

Polynomial: The polynomial is given as a (m+1) x (n+1) matrix/list of lists of (real) coefficients, in the example below the terms of the coefficients are given in their position:

[   1 * 1,   1 * x,   1 * x^2,   1 * x^3,  ... , 1 * x^n ]
[   y * 1,   y * x,   y * x^2,   y * x^4,  ... , y * x^n ]
[   ...  ,   ...   ,   ...   ,    ...   ,  ... ,   ...   ]
[ y^m * 1, y^m * x, y^m * x^2, y^m * x^3 , ..., y^m * x^n]


If you prefer, you can assume the matrix to be square (which can always be done with the necessary zero-padding), and if you want, you can also assume that the size of the matrix is given as aditional inputs.

In the following, the examples from above are represented as a matrix defined like this:

Circle:       Ellipse:      Parabola:  Cross:    Elliptic Curve: e.t.c
[-1, 0, 1]    [-1, 0, 1]    [ 0,-1]    [ 0, 0]   [-1, 1, 0,-1]
[ 0, 0, 0]    [ 0, 0, 0]    [ 0, 0]    [ 0, 1]   [ 0, 0, 0, 0]
[ 1, 0, 0]    [ 2, 0, 0]    [ 1, 0]              [ 1, 0, 0, 0]


Test cases with x-range / y-range:

(In a not so readable but better copy-paste-able format available here on pastebin.)

Circle:
[-1, 0, 1]   [-2,2]   [-2,2]
[ 0, 0, 0]
[ 1, 0, 0]

Ellipse:
[-1, 0, 1]   [-2,2]   [-1,1]
[ 0, 0, 0]
[ 2, 0, 0]

Cross:
[ 0, 0]      [-1,2]   [-2,1]
[ 0, 1]

Parabola:
[ 0,-1]      [-1,3]   [-2,2]
[ 0, 0]
[ 1, 0]

Elliptic Curve:
[-1, 1, 0,-1]    [-2,2]   [-3,3]
[ 0, 0, 0, 0]
[ 1, 0, 0, 0]

Folium of Descartes:
[  0,  0,  0,  1]    [-3,3]   [-3,3]
[  0, -3,  0,  0]
[  0,  0,  0,  0]
[  1,  0,  0,  0]

Lemniscate:
[  0,  0, -1,  0,  1]    [-2,2]   [-1,1]
[  0,  0,  0,  0,  0]
[  1,  0,  0,  0,  0]

Trifolium:
[ 0, 0, 0,-1, 1]    [-1,1]   [-1,1]
[ 0, 0, 0, 0, 0]
[ 0, 3, 2, 0, 0]
[ 0, 0, 0, 0, 0]
[ 1, 0, 0, 0, 0]

Astroid:
[ -1,  0,  3,  0, -3,  0,  1]    [-1,1]   [-1,1]
[  0,  0,  0,  0,  0,  0,  0]
[  3,  0, 21,  0,  3,  0,  0]
[  0,  0,  0,  0,  0,  0,  0]
[ -3,  0,  3,  0,  0,  0,  0]
[  0,  0,  0,  0,  0,  0,  0]
[  1,  0,  0,  0,  0,  0,  0]


I've got the inspiration for some curves from this pdf.

• Does "You do not have to worry about aliasing" mean that we can just colour each pixel according to whether its centre lies on the line? Jul 29, 2016 at 15:19
• I don't see the connection to aliasing. But no, there should be a continuous line separating the regions of different signs. Jul 29, 2016 at 15:32
• The matrix is not mxn, but (m+1)x(n+1). What do we take as input: m, n, or m+1,n+1? Or can we choose? Jul 29, 2016 at 23:18
• Can we just show the graphed function in a new window? Jul 30, 2016 at 7:35
• @LuisMendo Yes, the axis can be in any direction you like. (As long they are orthogonal=) Jul 30, 2016 at 18:13

The function g should be called with the matrix and the two ranges as arguments. The matrix is just a list of lists, the ranges each a two element list.

import Data.List
t=transpose
u=tail
z=zipWith
l%x=sum$z(*)l$iterate(*x)1                                   --generate powers and multiply with coefficients
e m y x=[l%x|l<-m]%y                                         --evaluate the encoded polynomial
a#b=[a,a+(b-a)/102..b]                                       --create a range
g m[u,w][i,j]=unlines$v[map((0<).e m y)$u#w|y<-i#j]          --evaluate the function on the grid, get the sign
f g=u[u$u$map fst$scanl(\(r,l)c->(c==l,c))(1<0,1<0) l|l<-g] --find +- or -+ transitions within lines a&b|a&&b=' '|0<1='#' --helper function for creating the string v g=z(z(&))(f g)(t$f$t g) --create the string  Here the outputs for the more interesting cases: Note that I had to downsize the resultion from 100x100 to about 40x40 such that it sill fits in the console (just change the hardcoded 102 to a smaller number). Also note that the y-axis is pointing downwards. • There are a couple of pretty small golfs you can make here. The last line uses parens when it could use $ to save a byte. Both places where you use map could be (<$>), and since you only use e once you can pull the (0<) inside it's definition. Also e could be named (!) to save 3 bytes. Dec 17, 2019 at 17:28 • And infixing z in the definition of v allows you to get rid of 4 parentheses (around z(&) and f g). Dec 17, 2019 at 17:33 • You could also rename # to a single character (e.g. s) and have it pattern match on the lists instead of g. (e.g. s[a,b]=[a,a+(b-a)/102..b];g m u i=unlines$v[m!y<\$>s u|y<-s i]) Dec 17, 2019 at 17:38

# Matlab, 114100 92 bytes

The right tool for the job? I use the interesting way Matlab does printf to generate a polynomial as a string. This polynomial can be provided to ezplot which plots the implicit curve on the specified domain. For readability the code is presented with newlines after ; which isn't needed and isn't counted towards the size.

function P(A,W,H,h,w)
t=0:h*w-1;
ezplot(sprintf('+%d*x^%.0f*y^%d',[A(:)';t/h;rem(t,h)]),[W,H])


## Golfing progress as expandable snippet.

Iteration 1, clean arguments, 114 bytes
---------------------------------------
function P(A,W,H)
h=size(A,1);
t=0:numel(A)-1;
ezplot(sprintf('+%d*x^%d*y^%d',[A(:)';floor(t/h);rem(t,h)]),[W,H])

Iteration 1.5, historical revisionism (thanks @flawr), 111 bytes
--------------------------------------------------------------
function P(A,W,H)
[h,w]=size(A);
t=0:h*w-1;
ezplot(sprintf('+%d*x^%d*y^%d',[A(:)';floor(t/h);rem(t,h)]),[W,H])

Iteration 2, playing dirty, 100 bytes
-------------------------------------
The code can be golfed further if I use that taking matrix sizes as inputs are
fair game, but it feels sad to ruin such neat input arguments. However,
this is code golf and I'll abuse anything I can and I end up at a nice 100 bytes:

function P(A,W,H,h,w)
t=0:h*w-1;
ezplot(sprintf('+%d*x^%d*y^%d',[A(:)';floor(t/h);rem(t,h)]),[W,H]);

Iteration 2.5, fix your code (thanks @LuisMendo), 98 bytes
----------------------------------------------------------
function P(A,W,H,h,w)
t=0:h*w-1;
ezplot(sprintf('+%d*x^%d*y^%d',[A(:)';fix(t/h);rem(t,h)]),[W,H]);

Iteration 3, more power to printf, 92 bytes [current best]
----------------------------------------------------------
Use %.0f instead of fix/floor.

function P(A,W,H,h,w)
t=0:h*w-1;
ezplot(sprintf('+%d*x^%.0f*y^%d',[A(:)';t/h;rem(t,h)]),[W,H])

Output of the test cases (click for full view): • Really nice solution using sprintf/ezplot! Jul 30, 2016 at 9:44
• Using fix instead of floor might help you reach the two-digit byte count :-) Jul 30, 2016 at 15:26
• You could also use [h,w]=size(A);t=0:h*w-1; to save another three bytes! Jul 30, 2016 at 18:16
• @LuisMendo Actually, I can do better. I was sad that Matlab's printf doesn't have a integer placeholder, but it still supports things like %.0f. That means I can drop flooring altogether and let printf fix it! Jul 30, 2016 at 18:19
• @flawr I do use the second part of that in later iterations. I realize my formatting with the best version last wasn't perfectly obvious. Edited formatting to make this crystal clear. Jul 30, 2016 at 18:41

# Python 2, 261 bytes

E=enumerate
M,[a,c],[b,d]=input()
e=(c-a)/199.
I=200
J=-int((b-d)/e-1)
print'P2',I,J,255
i=I*J
while i:i-=1;x,y=c-i%I*e,b+i/I*e;u,v,w=map(sum,zip(*((z*p/x,z*q/y,z)for q,R in E(M)for p,t in E(R)for z in[t*x**p*y**q])));print int(255*min(1,(w*w/(u*u+v*v))**.5/e))


Input format: matrix,xbounds,ybounds (e.g. [[-1,0,1],[0,0,0],[1,0,0]],[-2,2],[-2,2]). Output format: plain PGM.

This estimates the distance from every pixel center to the curve using the first-order approximation d(x, y) = |p(x, y)| / |∇p(x, y)|, where ∇p is the gradient of the polynomial p. (This is the distance from (x, y) to the intersection of the tangent plane at (x, y, p(x, y)) with the xy–plane.) Then the pixels where d(x, y) is less than one pixel width of the curve proportionally to d(x, y), resulting in nice antialiased lines (even though that isn’t a requirement). Here are the same graphs with the distance function divided by 16 to make it visible.

• Wait, so where in the code does the actual graphical plotting happen? Jul 30, 2016 at 7:50
• @R.Kap The code writes an image in plain PGM format to stdout. There’s one print statement for the image header and one print statement in the while loop for the value of each pixel. Jul 30, 2016 at 7:55
• Wow, that is really cool! Would you mind going a little bit more in depth about your plotting algorithm? Jul 30, 2016 at 9:48
• @flawr I’ve expanded the explanation a bit; does that answer your questions? Aug 2, 2016 at 0:20
• @AndersKaseorg Yes, thank you very much! Aug 2, 2016 at 9:53

# Python 3.5 + MatPlotLib + Numpy, 352 bytes:

from matplotlib.pyplot import*;from numpy import*
def R(M,S,U,r=range):N=linspace;E='+'.join([str(y)+'*'+m for y,m in[q for i,g in zip(M,[[i+'*'+p for p in['1']+['x^%d'%p for p in r(1,len(M))]]for i in['1']+['y^%d'%i for i in r(1,len(M))]])for q in zip(i,g)if q]]);x,y=meshgrid(N(*S,200),N(*U,200));contour(x,y,eval(E.replace('^','**')),0);show()


A named function. Pretty long, but hey, I'm just happy I was able to accomplish the task. Takes 3 inputs, which are the m by n matrix, the x-range, and y-range, which should all be in arrays (for example, [[-1,0,1],[0,0,0],[1,0,0]],[-2,2],[-2,2]). Outputs the completed graph in a new, graphical, interactive window. Will golf this down more of time when I can, but for now, I'm happy with it.

Final Outputs for the Test Cases: # MATL, 67 61 bytes

8Wt:qwq/t2:"wid*2M1)+i:q!^]!2&!w[1IK2]&!**ss&eZS5Y62&Y+|4=0YG


This code runs in release 18.5.0 of the language, which precedes the challenge. Input uses the optional m, n parameters. The matrix has semicolons as row separators. The exact input format (using the parabola as an example) is

[-1,3]
3
[-2,2]
2
[0,-1; 0, 0; 1, 0]


The code produces an image with size 255×255. This can be tested using @Suever's MATL Online compiler, which, among other very interesting features, includes graphical output. See for example

This compiler is still in an experimental stage. Please report any issues to @Suever in the MATL chatroom. If the "Run" button doesn't work, try refreshing the page and clicking again.

If you prefer ASCII output, the code needs to be modified a little (the changes only affect the first two and last four characters of the code above):

101t:qwq/t2:"wid*2M1)+i:q!^]!2&!w[1IK2]&!**ss&eZS5Y62&Y+|4<42*c


This produces a 100×100 ASCII grid which uses character * to represent the curve. You can also test this with @Dennis' Try it online! platform:

Note that the aspect ratio of ASCII output is altered because the characters are slightly higher than are wide.

### Explanation

The code first computes the two-variable polynomial on an x-y grid. This makes heavy use of broadcasting, computing an intermediate 4D array where each dimension represents x values, y values, x exponents, y exponents respectively.

From that function, the zero level line is computed. Since the challenge specifes that only sign changes need to be detected, the code applies 2D convolution with a 2×2 block of ones, and marks a pixel as belonging to the line if not the four values of the block have the same sign.

8W      % Push 2^8, that is, 256. (The ASCII-output version pushes 101 instead)
t:q     % Duplicate. Push range [0 1 ... 255]
wq      % Swap. Subtract 1 to obtain 255
/       % Divide. Gives normalized range [0 1/255 2/255... 1]
t       % Duplicate
2:"     % For loop: do this twice
w     %   Swap top two elements in the stack
i     %   Input two-number array defining x range (resp. y in second iteration)
d     %   Difference of the two entries
*     %   Multiply by normalized range
2M1)  %   Push the array again and get its first entry
+     %   Add. This gives the range for x values (resp. y)
i     %   Input m (n in second iteration)
:q    %   Range [0 1 ...m-1] (resp. [0 1 ...n-1])
!     %   Convert to column array
^     %   Power, element-wise with broadcast. This gives a matrix of size m×256
%   (resp. n×256) of powers of x (resp. y) for the range of values computed
%   previously
]       % End for loop
!       % Transpose. This transforms the n×256 matrix of powers of y into 256×n
2       % Push 2
&!      % Permute dimensions 1 and 3: transforms the 256×n matrix into a 4D array
% of size 1×n×256×1
w       % Swap top two elements in the stack: bring 256×m matrix to top
[1IK2]  % Push vector [1 3 4 2]
&!      % Permute dimensions as indicated by the vector: transforms the m×256 matrix
% into a 4D array of size m×1×1×256
*       % Multiply element-wise with broadcast: gives 4D array of size m×n×256×256
% with mixed powers of x and y for at the grid of x, y values
*       % Implicitly input m×n matrix. Multiply element-wise with broadcast: gives
% 4D array of size m×n×256×256
ss      % Sum along first two dimensions: gives 4D array of size 1×1×256×256
&e      % Squeeze singleton dimensions: gives matrix of size 256×256. This is the
% two-variable polynomial evaluated at the x, y grid.
% Now we need to find the zero level curve of this function. We do this by
% detecting when the sign of the function changes along any of the two axes
ZS      % Matrix of sign values (1, 0 or -1)
5Y6     % Predefined literal: matrix [1 1; 1 1]
2&Y+    % Compute 2D convolution, keeping only the valid (central) part
|4=     % True if absolute value of result is 4, which indicates no sign changes.
% (The ASCII version computes a negated version of this, for better display)
0YG     % Display as image. (The ASCII-output version does the following instead:
% multiply by 42 and convert to char. 42 is ASCII for '*', and character 0
% is shown as space. The 2D char array is then implicitly displayed)


### All test cases

Here are all inputs in the appropriate format, in case you want to try:

Circle:
[-2,2]
3
[-2,2]
3
[-1, 0, 1; 0, 0, 0; 1, 0, 0]

Ellipse:
[-2,2]
3
[-1,1]
3
[-1, 0, 1; 0, 0, 0; 2, 0, 0]

Cross:
[-1,2]
2
[-2,1]
2
[0, 0; 0, 1]

Parabola:
[-1,3]
3
[-2,2]
2
[0,-1; 0, 0; 1, 0]

Elliptic Curve:
[-2,2]
3
[-3,3]
4
[-1, 1, 0,-1; 0, 0, 0, 0; 1, 0, 0, 0]

Folium of Descartes:
[-3,3]
4
[-3,3]
4
[0,  0,  0,  1; 0, -3,  0,  0; 0,  0,  0,  0; 1,  0,  0,  0]

Lemniscate:
[-2,2]
3
[-1,1]
5
[0,  0, -1,  0,  1; 0,  0,  0,  0,  0; 1,  0,  0,  0,  0]

Trifolium:
[-1,1]
5
[-1,1]
5
[0, 0, 0,-1, 1; 0, 0, 0, 0, 0; 0, 3, 2, 0, 0; 0, 0, 0, 0, 0; 1, 0, 0, 0, 0]

Astroid
[-1,1]
7
[-1,1]
7
[-1,  0,  3,  0, -3,  0,  1; 0,  0,  0,  0,  0,  0,  0; 3,  0, 21,  0,  3,  0,  0; 0,  0,  0,  0,  0,  0,  0; -3,  0,  3,  0,  0,  0,  0; 0,  0,  0,  0,  0,  0,  0; 1,  0,  0,  0,  0,  0,  0]

• Still more readable than Perl. Great job, also a nice online compiler! Jul 30, 2016 at 19:04
• @flawr more readable than Perl LOL. As for the online compiler, it's all Suever's work! Jul 30, 2016 at 19:06
• @flawr Now with convolution! Jul 31, 2016 at 11:26
• <3 convolution! Jul 31, 2016 at 15:26