summaryrefslogtreecommitdiff
path: root/macros/luatex/latex/bezierplot
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /macros/luatex/latex/bezierplot
Initial commit
Diffstat (limited to 'macros/luatex/latex/bezierplot')
-rw-r--r--macros/luatex/latex/bezierplot/README21
-rw-r--r--macros/luatex/latex/bezierplot/bezierplot-doc.pdfbin0 -> 251276 bytes
-rw-r--r--macros/luatex/latex/bezierplot/bezierplot-doc.tex200
-rwxr-xr-xmacros/luatex/latex/bezierplot/bezierplot.lua940
-rw-r--r--macros/luatex/latex/bezierplot/bezierplot.sty17
5 files changed, 1178 insertions, 0 deletions
diff --git a/macros/luatex/latex/bezierplot/README b/macros/luatex/latex/bezierplot/README
new file mode 100644
index 0000000000..493bfe5c65
--- /dev/null
+++ b/macros/luatex/latex/bezierplot/README
@@ -0,0 +1,21 @@
+DESCRIPTION:
+bezierplot is a Lua program as well as a (Lua)LaTeX package.
+Given a smooth function, bezierplot returns a smooth bezier path written
+in TikZ notation (which also matches METAPOST) that approximates the
+graph of the function. For polynomial functions of degree <= 3 and
+inverses of them, the approximation is exact. bezierplot finds special
+points such as extreme points and inflection points and reduces the
+number of used points.
+
+VERSION:
+1.4 2019-06-18
+
+LICENSE:
+The package and the program are distributed on CTAN under the terms of
+the LaTeX Project Public License (LPPL) version 1.3c.
+
+Copyright (c) 2018 Linus Romer
+
+Please write to
+linus dot romer at gmx dot ch
+to submit bug reports, request new features, etc.
diff --git a/macros/luatex/latex/bezierplot/bezierplot-doc.pdf b/macros/luatex/latex/bezierplot/bezierplot-doc.pdf
new file mode 100644
index 0000000000..c93fd5f102
--- /dev/null
+++ b/macros/luatex/latex/bezierplot/bezierplot-doc.pdf
Binary files differ
diff --git a/macros/luatex/latex/bezierplot/bezierplot-doc.tex b/macros/luatex/latex/bezierplot/bezierplot-doc.tex
new file mode 100644
index 0000000000..a1239b6e57
--- /dev/null
+++ b/macros/luatex/latex/bezierplot/bezierplot-doc.tex
@@ -0,0 +1,200 @@
+% !TEX program = pdfLaTeX --shell-escape
+\documentclass[a4paper]{article}
+\usepackage{tikz,multicol,bezierplot,amsmath,cancel}
+\usepackage[margin=3.5cm,top=1.75cm]{geometry}
+\usepackage{fetamont}
+\title{bezierplot}\author{Linus Romer}
+\DeclareDocumentCommand{\graphcomparison}{ m m }{
+ \begin{center}
+ \begin{tikzpicture}[scale=.4]
+ \draw (0,-5) node[below]{\tiny\texttt{\detokenize{#1}\quad | \detokenize{#2}}};
+ \draw[step=1,thin] (-5,-5) grid (5,5);
+ \draw[thick,->] (-5,0) -- (5.5,0) node[below]{$x$};
+ \draw[thick,->] (0,-5) -- (0,5.5) node[left]{$y$};
+ \foreach \x in {-4,-3,-2,-1,1,2,3,4} {\draw (\x,1pt) -- (\x,-1pt) node[below]{\tiny \x};}
+ \foreach \y in {-4,-3,-2,-1,1,2,3,4} {\draw (1pt,\y) -- (-1pt,\y) node[left]{\tiny \y};}
+ \draw[color=red,domain=-5:5,range=-5:5,samples=1000] plot function{#2};
+ \draw \bezierplot{#1};
+ \end{tikzpicture}
+ \end{center}
+}
+\begin{document}
+\maketitle\noindent
+\section{Introduction}
+\texttt{bezierplot} is a Lua program as well as a (Lua)\LaTeX{} package. This document describes both.
+
+Given a smooth function, \texttt{bezierplot} returns a smooth bezier path written in Ti\emph{k}Z notation (which also matches \MP{}) that approximates the graph of the function. For polynomial functions of degree $\leq 3$ and inverses of them, the approximation is exact (up to numeric precision). \texttt{bezierplot} finds special graph points such as extreme points and inflection points and reduces the number of used points.
+
+The following example will show a comparison of \textsc{gnuplot} with \verb|bezierplot| for the function $y=\sqrt{x}$ for $0\leq x \leq 5$:
+\begin{center}
+ \begin{tikzpicture}[scale=1.4]
+ \draw (0,0) .. controls (0,0.7454) and (1.6667,1.4907) .. (5,2.2361);
+ \draw (0,0) circle(.02) -- (0,0.745) circle( .02);
+ \draw (1.6667,1.4907) circle(.02) -- (5,2.2361) circle( .02);
+ \draw (2.5,.5) node[above]{\verb|bezierplot|};
+ \begin{scope}[shift={(5.2,0)}]
+ \draw[domain=0:5,samples=51] plot function{x**0.5};
+ \foreach \x in {0,0.1,...,5.05} {\draw (\x,{\x^0.5}) circle (0.02);}
+ \draw (2.5,.5) node[above]{\textsc{gnuplot}};
+ \end{scope}
+ \end{tikzpicture}
+\end{center}
+\textsc{gnuplot} used 51 samples (no smoothing) and is still quite inexact at the beginning, whereas \verb|bezierplot| uses 4 points only and is exact (up to numeric precision)!
+\section{Installation}
+As \texttt{bezierplot} is written in Lua, the installation depends whether you are using Lua\LaTeX{} or another \LaTeX{} engine.
+\subsection{Installation For Lua\LaTeX{}}
+If you have installed \texttt{bezierplot} by a package manager, the installation is already complete. The manual installation of \texttt{bezierplot} is done in 2 steps:
+\begin{itemize}
+ \item copy the files \texttt{bezierplot.lua} and \texttt{bezierplot.sty} somewhere in your \texttt{texmf} tree (e.g. to \verb|~/texmf/tex/lualatex/bezierplot/bezierplot.sty| and\\
+ \verb|~/texmf/scripts/bezierplot/bezierplot.lua|)
+ \item update the ls-R databases by running \texttt{mktexlsr}
+\end{itemize}
+\subsection{Additional Installation Steps For Other \LaTeX{} Engines}
+You will have to call \texttt{bezierplot} as an external program via the option \texttt{--shell-escape} (\texttt{--write18} for MiK\TeX{}). Therefore, \texttt{bezierplot.lua} has to be copied with the name \texttt{bezierplot} to a place, where your OS can find it. Under Linux this usually means copying to the directory \texttt{/usr/local/bin/}, but for Windows this will probably include more steps (like adding to the \texttt{PATH}). Of course, Lua has to be installed as well. As soon as you can call \texttt{bezierplot} from a command line (e.g. by typing \verb|bezierplot "x^2"|), it should also work with other \LaTeX{} engines.
+\section{Loading}
+The \texttt{bezierplot} package is loaded with \verb|\usepackage{bezierplot}|. There are no loading options for the package.
+\section{Usage}
+\begin{multicols}{2}
+\noindent A minimal example of Lua\LaTeX{} document could be:
+\begin{verbatim}
+\documentclass{article}
+\usepackage{tikz,bezierplot}
+\begin{document}
+\tikz \draw \bezierplot{x^2};
+\end{document}
+\end{verbatim}
+\begin{center}
+ \tikz \draw[scale=.7] \bezierplot{x^2};
+\end{center}
+\end{multicols}
+\noindent
+The command \verb|\bezierplot| has 6 optional arguments in the sense of
+\begin{center}
+ \verb|\bezierplot[XMIN][XMAX][YMIN][YMAX][SAMPLES]{FUNCTION}|
+\end{center}
+The defaults are \verb|XMIN| = \verb|YMIN| $= -5$, \verb|XMAX| = \verb|YMAX| $= 5$ and \verb|SAMPLES| $= 0$ (this will set as few samples as possible).
+\begin{center}
+ \begin{tikzpicture}[scale=.7]
+ \draw \bezierplot[-1][2]{x^2};
+ \draw (0,0) node[below]{\verb|\bezierplot[-1][2]{x^2}|};
+ \begin{scope}[shift={(10,0)}]
+ \draw \bezierplot[-1][2][0.5][3]{x^2};
+ \draw (0,0) node[below]{\verb|\bezierplot[-1][2][0.5][3]{x^2}|};
+ \end{scope}
+ \end{tikzpicture}
+\end{center}
+You may reverse the graph by making \verb|XMIN| bigger than \verb|XMAX|. E.g.
+\begin{verbatim}
+ \bezierplot[-5][5]{0.5*x+1}
+\end{verbatim}
+returns \verb|(-5,-1.5) -- (5,3.5)|, whereas
+\begin{verbatim}
+ \bezierplot[5][-5]{0.5*x+1}
+\end{verbatim}
+returns the reversed path \verb|(5,3.5) -- (-5,-1.5)|. This is useful, if you want to cycle a path to a closed area:
+\begin{multicols}{2}
+\begin{verbatim}
+\begin{tikzpicture}
+ \fill[black!30] \bezierplot[-1][1]{2-x^2}
+ -- \bezierplot[1][-1]{x^3-x} -- cycle;
+ \draw \bezierplot[-1.1][1.1]{2-x^2};
+ \draw \bezierplot[-1.1][1.1]{x^3-x};
+\end{tikzpicture}
+\end{verbatim}
+\begin{center}
+ \begin{tikzpicture}
+ \fill[black!30] \bezierplot[-1][1]{2-x^2} -- \bezierplot[1][-1]{x^3-x} -- cycle;
+ \draw \bezierplot[-1.1][1.1]{2-x^2};
+ \draw \bezierplot[-1.1][1.1]{x^3-x};
+ \end{tikzpicture}
+\end{center}
+\end{multicols}
+\subsection{Running Raw \texttt{bezierplot}}
+Of course, you can run \verb|bezierplot.lua| in a terminal without using \LaTeX{}, e.g.
+\begin{verbatim}
+lua bezierplot.lua "3*x^0.8+2"
+\end{verbatim}
+will return
+\begin{verbatim}
+(0,2) .. controls (0.03,2.282) and (0.268,3.244) .. (1,5)
+\end{verbatim}
+You can set the window of the graph and the number of samples as follows:
+\begin{verbatim}
+lua bezierplot.lua "FUNCTION" XMIN XMAX YMIN YMAX SAMPLES
+\end{verbatim}
+e.g.
+\begin{verbatim}
+lua bezierplot.lua "FUNCTION" 0 1 -3 2.5 201
+\end{verbatim}
+will set $0\leq x\leq 1$ and $-3\leq y\leq 2.5$ and $201$ equidistant samples. You may also omit the $y$--range, hence
+\begin{verbatim}
+lua bezierplot.lua "FUNCTION" 0 1
+\end{verbatim}
+will set $0\leq x\leq 1$ and leave the default $-5\leq y\leq 5$. The variables \verb|XMIN|, \verb|XMAX|, \verb|YMIN| and \verb|YMAX| may also be computable expressions like \verb|2*pi+6|:
+\begin{verbatim}
+lua bezierplot.lua "sin(x)" -pi pi
+\end{verbatim}
+You may use \verb|huge| for $\infty$:
+\begin{verbatim}
+lua bezierplot "1/x" 0 1 0 huge
+\end{verbatim}
+As \verb|huge| is very huge and \verb|bezierplot| uses recursive calls for nontrivial functions and non--fixed samples, this can last very long:
+\begin{verbatim}
+lua bezierplot "1/x" -5 5 -huge huge
+\end{verbatim}
+But if you set fixed samples, it will be fast again (as this does not use recursive calls):
+\begin{verbatim}
+lua bezierplot "1/x" -5 5 -huge huge 100
+\end{verbatim}
+\subsection{Notation Of Functions}
+The function term given to \verb|bezierplot| must contain at most one variable: $x$. E.g. \verb|"2.3*(x-1)^2-3"|. You must not omit \verb|*| operators:
+\begin{center}
+ wrong:\quad $\cancel{\texttt{2x(x+1)}}$ \hfil correct:\quad \texttt{2*x*(x+1)}
+\end{center}
+You have two possibilities to write powers: \verb|"x^2"| and \verb|"x**2"| both mean $x^2$.
+
+\medskip
+
+The following functions and constants are possible:
+\begin{center}
+\begin{tabular}{ll}
+ \verb|abs| & absolute value (remember: your function should still be smooth)\\
+ \verb|acos| & $\cos^{-1}$ inverse function of cosine in radians\\
+ \verb|asin| & $\sin^{-1}$ inverse function of sine in radians\\
+ \verb|atan| & $\tan^{-1}$ inverse function of tangent in radians\\
+ \verb|cbrt| & cube root $\sqrt[3]{\quad}$ that works for negative numbers, too\\
+ \verb|cos| & cosine for angles in radians\\
+ \verb|exp| & the exponential function $e^{(\;)}$\\
+ \verb|huge| & the numerical $\infty$\\
+ \verb|e| & the euler constant $e=\mathrm{exp}(1)$\\
+ \verb|log| & the natural logarithm $\mathrm{log}_e(\;)$\\
+ \verb|pi| & Archimedes’ constant $\pi\approx 3.14$\\
+ \verb|sgn| & sign function\\
+ \verb|sin| & sine for angles in radians\\
+ \verb|sqrt| & square root $\sqrt{\quad}$\\
+ \verb|tan| & tangent for angles in radians
+\end{tabular}
+\end{center}
+%
+\newpage
+%
+\section{Examples of \texttt{bezierplot} in Comparison with \textsc{gnuplot}}
+The following graphs are drawn with \texttt{bezierplot} (black) and \textsc{gnuplot} (red). You may not recognize the red behind the black unless you zoom in. \textsc{gnuplot} used 1000 samples per example. The functions are given below the pictures (left: bezierplot, right: \textsc{gnuplot}).
+\begin{multicols}{3}
+\graphcomparison{0.32*x-0.7}{0.32*x-0.7}
+\graphcomparison{-x^2+4}{-x**2+4}
+\graphcomparison{(x+1)*x*(x-1)}{(x+1)*x*(x-1)}
+\graphcomparison{x^0.5}{x**0.5}
+%\graphcomparison{x^(1/3)}{x**(1/3.)}
+\graphcomparison{cbrt(x)}{sgn(x)*abs(x)**(1/3.)}
+\graphcomparison{x^3*(x-1)}{x**3*(x-1)}
+\graphcomparison{2*cos(3*x+4)+3}{2*cos(3*x+4)+3}
+\graphcomparison{tan(x)}{tan(x)}
+\graphcomparison{x+0.5*sin(x)}{x+0.5*sin(x)}
+%\graphcomparison{1/(x-2)+1}{1/(x-2)+1}
+\graphcomparison{2*x^2/(3*x-3)}{2*x**2/(3*x-3)}
+\graphcomparison{4-exp(x)}{4-exp(x)}
+\graphcomparison{log(x+4)}{log(x+4)}
+\end{multicols}
+
+\end{document}
diff --git a/macros/luatex/latex/bezierplot/bezierplot.lua b/macros/luatex/latex/bezierplot/bezierplot.lua
new file mode 100755
index 0000000000..382da2d74e
--- /dev/null
+++ b/macros/luatex/latex/bezierplot/bezierplot.lua
@@ -0,0 +1,940 @@
+#!/usr/bin/env lua
+-- Linus Romer, published 2018 under LPPL Version 1.3c
+-- version 1.4 2019-06-18
+abs = math.abs
+acos = math.acos
+asin = math.asin
+atan = math.atan
+cos = math.cos
+exp = math.exp
+e = math.exp(1)
+log = math.log
+pi = math.pi
+sin = math.sin
+sqrt = math.sqrt
+tan = math.tan
+huge = math.huge
+
+-- just a helper for debugging:
+local function printdifftable(t)
+ for i = 1,#t do
+ for j = 1, 6 do
+ if j < 5 then
+ io.write(t[i][j].." ")
+ else
+ if t[i][j] then
+ io.write("true ")
+ else
+ io.write("false ")
+ end
+ end
+ end
+ io.write("\n")
+ end
+end
+
+-- cube root defined for all real numbers x
+function cbrt(x)
+ if x < 0 then
+ return -(-x)^(1/3)
+ else
+ return x^(1/3)
+ end
+end
+
+function sgn(x)
+ if x<0 then
+ return -1
+ elseif x>0 then
+ return 1
+ else
+ return 0
+ end
+end
+
+local function round(num, decimals)
+ local result = tonumber(string.format("%." .. (decimals or 0) .. "f", num))
+ if abs(result) == 0 then
+ return 0
+ else
+ return result
+ end
+end
+
+-- check if a point (x,y) satisfies xmin <= x <= xmax and ymin <= < <= ymax
+local function is_in_window(x,y,xmin,xmax,ymin,ymax)
+ if x >= xmin and x <= xmax and y >= ymin and y <= ymax then
+ return true
+ else
+ return false
+ end
+end
+
+local function evaluate(s)
+ local tempfunc = assert(load("return " .. s))
+ return tempfunc()
+end
+
+
+-- 5-stencil method
+-- return from a graph from f in the form {{x,y},...}
+-- the derivatives in form {{x,y,dy/dx,ddy/ddx,extrema,inflection},...}
+local function diffgraph(func,graph,h)
+ local dgraph = {}
+ local l = #graph
+ if l < 4 then -- this is not worth the pain...
+ for i = 1, l do
+ table.insert(dgraph,{graph[i][1],graph[i][2],0,0,false,false})
+ end
+ else
+ local yh = func(graph[1][1]-h)
+ local yhh = func(graph[1][1]-2*h)
+ if yhh > -math.huge and yhh < math.huge -- if defined at all
+ and yh > -math.huge and yh < math.huge then
+ dgraph[1] = {graph[1][1],graph[1][2],
+ (yhh-8*yh+8*graph[2][2]-graph[3][2])/(12*h),
+ (-yhh+16*yh-30*graph[1][2]+16*graph[2][2]-graph[3][2])
+ /(12*h^2)}
+ dgraph[2] = {graph[2][1],graph[2][2],
+ (yh-8*graph[1][2]+8*graph[3][2]-graph[4][2])/(12*h),
+ (-yh+16*graph[1][2]-30*graph[2][2]+16*graph[3][2]-graph[4][2])
+ /(12*h^2)}
+ else -- take neighbour values
+ dgraph[1] = {graph[1][1],graph[1][2],
+ (graph[1][2]-8*graph[2][2]+8*graph[4][2]-graph[5][2])/(12*h),
+ (-graph[1][2]+16*graph[2][2]-30*graph[3][2]
+ +16*graph[4][2]-graph[5][2])/(12*h^2)}
+ dgraph[2] = {graph[2][1],graph[2][2],
+ (graph[1][2]-8*graph[2][2]+8*graph[4][2]-graph[5][2])/(12*h),
+ (-graph[1][2]+16*graph[2][2]-30*graph[3][2]
+ +16*graph[4][2]-graph[5][2])/(12*h^2)}
+ end
+ for i = 3, l-2 do
+ table.insert(dgraph,{graph[i][1],graph[i][2],
+ (graph[i-2][2]-8*graph[i-1][2]+8*graph[i+1][2]-graph[i+2][2])
+ /(12*h),
+ (-graph[i-2][2]+16*graph[i-1][2]-30*graph[i][2]
+ +16*graph[i+1][2]-graph[i+2][2])
+ /(12*h^2)})
+ end
+ yh = func(graph[l][1]+h)
+ yhh = func(graph[l][1]+2*h)
+ if yhh > -math.huge and yhh < math.huge -- if defined at all
+ and yh > -math.huge and yh < math.huge then
+ dgraph[l-1] = {graph[l-1][1],graph[l-1][2],
+ (graph[l-3][2]-8*graph[l-2][2]+8*graph[l][2]-yh)/(12*h),
+ (-graph[l-3][2]+16*graph[l-2][2]-30*graph[l-1][2]
+ +16*graph[l][2]-yh)/(12*h^2)}
+ dgraph[l] = {graph[l][1],graph[l][2],
+ (graph[l-2][2]-8*graph[l-1][2]+8*yh-yhh)/(12*h),
+ (-graph[l-2][2]+16*graph[l-1][2]-30*graph[l][2]
+ +16*yh-yhh)/(12*h^2)}
+ else
+ -- take neighbour values
+ dgraph[l] = {graph[l][1],graph[l][2],
+ (graph[l-4][2]-8*graph[l-3][2]+8*graph[l-1][2]-graph[l][2])
+ /(12*h),
+ (-graph[l-4][2]+16*graph[l-3][2]-30*graph[l-2][2]
+ +16*graph[l-1][2]-graph[l][2])/(12*h^2)}
+ dgraph[l-1] = {graph[l-1][1],graph[l-2][2],
+ (graph[l-4][2]-8*graph[l-3][2]+8*graph[l-1][2]-graph[l][2])
+ /(12*h),
+ (-graph[l-4][2]+16*graph[l-3][2]-30*graph[l-2][2]
+ +16*graph[l-1][2]-graph[l][2])/(12*h^2)}
+ end
+ -- add information about being extremum / inflection point (true/false)
+ for i = 1, l do
+ dgraph[i][5] = false -- dy/dx == 0 ? default, may change later
+ dgraph[i][6] = false -- ddy/ddx == 0 ? default, may change later
+ end
+ for i = 1, l-1 do
+ -- if no gap is inbetween
+ if not (dgraph[i+1][1] - dgraph[i][1] > 1.5*h) then
+ -- check for dy/dx == 0
+ -- if not already determined as near dy/dx=0
+ if dgraph[i][3] == 0 then
+ if dgraph[i+1][3] == 0 then --take the later
+ dgraph[i+1][5] = true
+ dgraph[i][5] = false
+ else
+ dgraph[i][5] = true
+ end
+ elseif abs(dgraph[i][3]*dgraph[i+1][3])
+ ~= dgraph[i][3]*dgraph[i+1][3] then -- this must be near
+ if abs(dgraph[i][4]) <= abs(dgraph[i+1][4]) then
+ dgraph[i][5] = true
+ else
+ dgraph[i+1][5] = true
+ end
+ end
+ -- check for ddy/ddx == 0
+ -- if not already determined as near ddy/ddx=0
+ if (not dgraph[i][6]) and (abs(dgraph[i][4]*dgraph[i+1][4])
+ ~= dgraph[i][4]*dgraph[i+1][4]) then -- this must be near
+ if abs(dgraph[i][4]) <= abs(dgraph[i+1][4]) then
+ dgraph[i][6] = true
+ else
+ dgraph[i+1][6] = true
+ end
+ end
+ end
+ end
+ end
+ return dgraph
+end
+
+-- simplified diffgraph function, the function is derived only once
+-- return from a graph from f in the form {{x,y},...}
+-- the derivatives in form {{x,y,dy/dx},...}
+-- we start with index 1 and then always jump indexjump to the next
+-- index
+local function diffgraphsimple(func,graph,h,indexjump)
+ local dgraph = {}
+ local l = math.floor((#graph-1)/indexjump)*indexjump+1
+ if l < 2 then -- this is not worth the pain...
+ dgraph = {graph[1][1],graph[1][2],0}
+ else
+ local yh = func(graph[1][1]-h)
+ local yhh = func(graph[1][1]-2*h)
+ if yhh > -math.huge and yhh < math.huge -- if defined at all
+ and yh > -math.huge and yh < math.huge then
+ dgraph[1] = {graph[1][1],graph[1][2],
+ (yhh-8*yh+8*graph[2][2]-graph[3][2])/(12*h)}
+ else -- take neighbour values
+ dgraph[1] = {graph[1][1],graph[1][2],
+ (graph[1][2]-8*graph[2][2]+8*graph[4][2]-graph[5][2])
+ /(12*h)}
+ end
+ for i = 1+indexjump, l-1, indexjump do
+ table.insert(dgraph,{graph[i][1],graph[i][2],
+ (graph[i-2][2]-8*graph[i-1][2]+8*graph[i+1][2]-graph[i+2][2])
+ /(12*h)})
+ end
+ yh = func(graph[l][1]+h)
+ yhh = func(graph[l][1]+2*h)
+ if yhh > -math.huge and yhh < math.huge -- if defined at all
+ and yh > -math.huge and yh < math.huge then
+
+ table.insert(dgraph,{graph[l][1],graph[l][2],
+ (graph[l-2][2]-8*graph[l-1][2]+8*yh-yhh)/(12*h)})
+ else
+ -- take neighbour values
+ table.insert(dgraph,{graph[l][1],graph[l][2],
+ (graph[l-4][2]-8*graph[l-3][2]+8*graph[l-1][2]-graph[l][2])
+ /(12*h)})
+ end
+ end
+ return dgraph
+end
+
+-- diffgraph for cubic function ax^3+bx^2+cx+d
+-- returns the derivatives in form {{x,y,dy/dx,ddy/ddx},...}
+-- if isinverse = true then the coordinates will be inversed
+local function diffgraphcubic(graph,a,b,c,d,isinverse)
+ local dgraph = {}
+ local l = #graph
+ for i = 1, l do
+ if isinverse then
+ dgraph[#dgraph+1] = {graph[i][2],graph[i][1],c
+ +graph[i][2]*(2*b+3*a*graph[i][2]),6*a*graph[i][2]+2*b}
+ else
+ dgraph[#dgraph+1] = {graph[i][1],graph[i][2],c
+ +graph[i][1]*(2*b+3*a*graph[i][1]),6*a*graph[i][1]+2*b}
+ end
+ end
+ return dgraph
+end
+
+-- checks for 100 x, if the function given by funcstring
+-- fits the graph g (up to maxerror) after filling in
+-- the parameters a, b, c, d
+-- if the graph is inverted, then isinverse has to be set true
+local function do_parameters_fit(a,b,c,d,funcstring,funcgraph,maxerror,isinverse)
+ if not (a > -math.huge and a < math.huge and b > -math.huge and b < math.huge and
+ c > -math.huge and c < math.huge and d > -math.huge and d < math.huge) then
+ return false
+ end
+ local funcx = string.gsub(funcstring, "a", a)
+ local funcx = string.gsub(funcx, "b", b)
+ local funcx = string.gsub(funcx, "c", c)
+ local funcx = string.gsub(funcx, "d", d)
+ local func = assert(load("local x = ...; return "..funcx))
+ for i = 1, #funcgraph, math.max(1,math.floor(0.01*#funcgraph)) do
+ if isinverse then
+ if abs(func(funcgraph[i][2])-funcgraph[i][1])
+ > maxerror then
+ return false
+ end
+ else
+ if abs(func(funcgraph[i][1])-funcgraph[i][2])
+ > maxerror then
+ return false
+ end
+ end
+ end
+ return true
+end
+
+-- f(x)=a*x^3+b*x+c
+local function parameters_cubic(xp,yp,xq,yq,xr,yr,xs,ys)
+ local a = (((xp^2 * xq) * yr) - ((xp^2 * xq) * ys)
+ - ((xp^2 * xr) * yq) + ((xp^2 * xr) * ys) + ((xp^2 * xs) * yq)
+ - ((xp^2 * xs) * yr) - ((xp * xq^2) * yr) + ((xp * xq^2) * ys)
+ + ((xp * xr^2) * yq) - ((xp * xr^2) * ys) - ((xp * xs^2) * yq)
+ + ((xp * xs^2) * yr) + ((xq^2 * xr) * yp) - ((xq^2 * xr) * ys)
+ - ((xq^2 * xs) * yp) + ((xq^2 * xs) * yr) - ((xq * xr^2) * yp)
+ + ((xq * xr^2) * ys) + ((xq * xs^2) * yp) - ((xq * xs^2) * yr)
+ + ((xr^2 * xs) * yp) - ((xr^2 * xs) * yq) - ((xr * xs^2) * yp)
+ + ((xr * xs^2) * yq)) /
+ (((xp^3 * xq^2) * xr) - ((xp^3 * xq^2) * xs)
+ - ((xp^3 * xq) * xr^2) + ((xp^3 * xq) * xs^2)
+ + ((xp^3 * xr^2) * xs) - ((xp^3 * xr) * xs^2)
+ - ((xp^2 * xq^3) * xr) + ((xp^2 * xq^3) * xs)
+ + ((xp^2 * xq) * xr^3) - ((xp^2 * xq) * xs^3)
+ - ((xp^2 * xr^3) * xs) + ((xp^2 * xr) * xs^3)
+ + ((xp * xq^3) * xr^2) - ((xp * xq^3) * xs^2)
+ - ((xp * xq^2) * xr^3) + ((xp * xq^2) * xs^3)
+ + ((xp * xr^3) * xs^2) - ((xp * xr^2) * xs^3)
+ - ((xq^3 * xr^2) * xs) + ((xq^3 * xr) * xs^2)
+ + ((xq^2 * xr^3) * xs) - ((xq^2 * xr) * xs^3)
+ - ((xq * xr^3) * xs^2) + ((xq * xr^2) * xs^3))
+ local b = ((((-xp^3) * xq) * yr) + ((xp^3 * xq) * ys)
+ + ((xp^3 * xr) * yq) - ((xp^3 * xr) * ys) - ((xp^3 * xs) * yq)
+ + ((xp^3 * xs) * yr) + ((xp * xq^3) * yr) - ((xp * xq^3) * ys)
+ - ((xp * xr^3) * yq) + ((xp * xr^3) * ys) + ((xp * xs^3) * yq)
+ - ((xp * xs^3) * yr) - ((xq^3 * xr) * yp) + ((xq^3 * xr) * ys)
+ + ((xq^3 * xs) * yp) - ((xq^3 * xs) * yr) + ((xq * xr^3) * yp)
+ - ((xq * xr^3) * ys) - ((xq * xs^3) * yp) + ((xq * xs^3) * yr)
+ - ((xr^3 * xs) * yp) + ((xr^3 * xs) * yq) + ((xr * xs^3) * yp)
+ - ((xr * xs^3) * yq)) /
+ (((xp^3 * xq^2) * xr) - ((xp^3 * xq^2) * xs)
+ - ((xp^3 * xq) * xr^2) + ((xp^3 * xq) * xs^2)
+ + ((xp^3 * xr^2) * xs) - ((xp^3 * xr) * xs^2)
+ - ((xp^2 * xq^3) * xr) + ((xp^2 * xq^3) * xs)
+ + ((xp^2 * xq) * xr^3) - ((xp^2 * xq) * xs^3)
+ - ((xp^2 * xr^3) * xs) + ((xp^2 * xr) * xs^3)
+ + ((xp * xq^3) * xr^2) - ((xp * xq^3) * xs^2)
+ - ((xp * xq^2) * xr^3) + ((xp * xq^2) * xs^3)
+ + ((xp * xr^3) * xs^2) - ((xp * xr^2) * xs^3)
+ - ((xq^3 * xr^2) * xs) + ((xq^3 * xr) * xs^2)
+ + ((xq^2 * xr^3) * xs) - ((xq^2 * xr) * xs^3)
+ - ((xq * xr^3) * xs^2) + ((xq * xr^2) * xs^3))
+ local c = (((xp^3 * xq^2) * yr) - ((xp^3 * xq^2) * ys)
+ - ((xp^3 * xr^2) * yq) + ((xp^3 * xr^2) * ys)
+ + ((xp^3 * xs^2) * yq) - ((xp^3 * xs^2) * yr)
+ - ((xp^2 * xq^3) * yr) + ((xp^2 * xq^3) * ys)
+ + ((xp^2 * xr^3) * yq) - ((xp^2 * xr^3) * ys)
+ - ((xp^2 * xs^3) * yq) + ((xp^2 * xs^3) * yr)
+ + ((xq^3 * xr^2) * yp) - ((xq^3 * xr^2) * ys)
+ - ((xq^3 * xs^2) * yp) + ((xq^3 * xs^2) * yr)
+ - ((xq^2 * xr^3) * yp) + ((xq^2 * xr^3) * ys)
+ + ((xq^2 * xs^3) * yp) - ((xq^2 * xs^3) * yr)
+ + ((xr^3 * xs^2) * yp) - ((xr^3 * xs^2) * yq)
+ - ((xr^2 * xs^3) * yp) + ((xr^2 * xs^3) * yq)) /
+ (((xp^3 * xq^2) * xr) - ((xp^3 * xq^2) * xs)
+ - ((xp^3 * xq) * xr^2) + ((xp^3 * xq) * xs^2)
+ + ((xp^3 * xr^2) * xs) - ((xp^3 * xr) * xs^2)
+ - ((xp^2 * xq^3) * xr) + ((xp^2 * xq^3) * xs)
+ + ((xp^2 * xq) * xr^3) - ((xp^2 * xq) * xs^3)
+ - ((xp^2 * xr^3) * xs) + ((xp^2 * xr) * xs^3)
+ + ((xp * xq^3) * xr^2) - ((xp * xq^3) * xs^2)
+ - ((xp * xq^2) * xr^3) + ((xp * xq^2) * xs^3)
+ + ((xp * xr^3) * xs^2) - ((xp * xr^2) * xs^3)
+ - ((xq^3 * xr^2) * xs) + ((xq^3 * xr) * xs^2)
+ + ((xq^2 * xr^3) * xs) - ((xq^2 * xr) * xs^3)
+ - ((xq * xr^3) * xs^2) + ((xq * xr^2) * xs^3))
+ local d = ((((xp^(3) * xq^(2)) * xr) * ys)
+ - (((xp^(3) * xq^(2)) * xs) * yr) - (((xp^(3) * xq) * xr^(2)) * ys)
+ + (((xp^(3) * xq) * xs^(2)) * yr) + (((xp^(3) * xr^(2)) * xs) * yq)
+ - (((xp^(3) * xr) * xs^(2)) * yq) - (((xp^(2) * xq^(3)) * xr) * ys)
+ + (((xp^(2) * xq^(3)) * xs) * yr) + (((xp^(2) * xq) * xr^(3)) * ys)
+ - (((xp^(2) * xq) * xs^(3)) * yr) - (((xp^(2) * xr^(3)) * xs) * yq)
+ + (((xp^(2) * xr) * xs^(3)) * yq) + (((xp * xq^(3)) * xr^(2)) * ys)
+ - (((xp * xq^(3)) * xs^(2)) * yr) - (((xp * xq^(2)) * xr^(3)) * ys)
+ + (((xp * xq^(2)) * xs^(3)) * yr) + (((xp * xr^(3)) * xs^(2)) * yq)
+ - (((xp * xr^(2)) * xs^(3)) * yq) - (((xq^(3) * xr^(2)) * xs) * yp)
+ + (((xq^(3) * xr) * xs^(2)) * yp) + (((xq^(2) * xr^(3)) * xs) * yp)
+ - (((xq^(2) * xr) * xs^(3)) * yp) - (((xq * xr^(3)) * xs^(2)) * yp)
+ + (((xq * xr^(2)) * xs^(3)) * yp)) /
+ (((xp^(3) * xq^(2)) * xr) -
+ ((xp^(3) * xq^(2)) * xs) - ((xp^(3) * xq) * xr^(2))
+ + ((xp^(3) * xq) * xs^(2)) + ((xp^(3) * xr^(2)) * xs)
+ - ((xp^(3) * xr) * xs^(2)) - ((xp^(2) * xq^(3)) * xr)
+ + ((xp^(2) * xq^(3)) * xs) + ((xp^(2) * xq) * xr^(3))
+ - ((xp^(2) * xq) * xs^(3)) - ((xp^(2) * xr^(3)) * xs)
+ + ((xp^(2) * xr) * xs^(3)) + ((xp * xq^(3)) * xr^(2))
+ - ((xp * xq^(3)) * xs^(2)) - ((xp * xq^(2)) * xr^(3))
+ + ((xp * xq^(2)) * xs^(3)) + ((xp * xr^(3)) * xs^(2))
+ - ((xp * xr^(2)) * xs^(3)) - ((xq^(3) * xr^(2)) * xs)
+ + ((xq^(3) * xr) * xs^(2)) + ((xq^(2) * xr^(3)) * xs)
+ - ((xq^(2) * xr) * xs^(3)) - ((xq * xr^(3)) * xs^(2))
+ + ((xq * xr^(2)) * xs^(3)))
+ return a, b, c, d
+end
+
+-- f(x)=a*x+b
+local function parameters_affine(xp,yp,xq,yq)
+ local a = (yp - yq) / (xp - xq)
+ local b = ((xp * yq) - (xq * yp)) / (xp - xq)
+ return a, b
+end
+
+-- what is the sum of the squared error
+-- when comparing the bezier path
+-- p.. control q and r .. s
+-- with the graph g from index starti to endi
+-- (looking at the points at roughly t=.33 and t=.67)
+local function squareerror(f,g,starti,endi,qx,qy,rx,ry)
+ local result = 0
+ for t = .1, .9, .1 do
+ x = (1-t)^3*g[starti][1]+3*t*(1-t)^2*qx+3*t^2*(1-t)*rx+t^3*g[endi][1]
+ y = (1-t)^3*g[starti][2]+3*t*(1-t)^2*qy+3*t^2*(1-t)*ry+t^3*g[endi][2]
+ result = result + (y-f(x))^2
+ end
+ return result
+end
+
+-- converts a table with bezier point information
+-- to a string with rounded values
+-- the path is reversed, if rev is true
+-- e.g. if bezierpoints = {{0,1},{2,3,4,5,6,7},{8,9,10,11,12,13}}
+-- then
+-- (0,1) .. controls (2,3) and (4,5) .. (6,7) .. controls
+-- (8,9) and (10,11) .. (12,13)
+-- will be returned
+-- the notation "pgfplots" will change the notation to
+-- YES: \addplot coordinates {(0,1) (6,7) (2,3) (4,5) (6,7) (12,13) (8,9) (10,11)}
+-- NO: 0 1 \\ 6 7 \\ 2 3 \\ 4 5 \\ \\ 6 7 \\ 12 13 \\ 8 9 \\ 10 11 \\
+-- As pgfplots does not connect the bezier segments
+-- reverse paths are not implemented
+local function beziertabletostring(bezierpoints,rndx,rndy,rev,notation)
+ local bezierstring = ""
+ local b = {{round(bezierpoints[1][1],rndx),round(bezierpoints[1][2],rndy)}} -- rounded and then
+ -- reduced points (if identical after rounding)
+ -- rounding
+ for i = 2, #bezierpoints do
+ -- check if x--coordinates are identical
+ if round(bezierpoints[i][#bezierpoints[i]-1],rndx) ~= b[#b][#b[#b]-1] then
+ b[#b+1] = {}
+ for j = 1, #bezierpoints[i] do
+ if j % 2 == 0 then -- x coordinate
+ b[#b][j] = round(bezierpoints[i][j],rndx)
+ else
+ b[#b][j] = round(bezierpoints[i][j],rndy)
+ end
+ end
+ end
+ end
+ if #b > 1 then -- if not empty or single point
+ -- check if bezierstring contains only straight lines
+ local onlystraightlines = true
+ for i = 1, #b do
+ if #b[i] > 2 then
+ onlystraightlines = false
+ break
+ end
+ end
+ if onlystraightlines then
+ if rev then
+ bezierstring = "(" .. b[#b][1] .. "," .. b[#b][2] ..")"
+ for i = #b-1, 1, -1 do
+ bezierstring = bezierstring ..
+ " -- (" .. b[i][1] .. "," .. b[i][2] ..")"
+ end
+ else
+ if notation == "pgfplots" then
+ bezierstring = "\\addplot coordinates {("
+ .. b[1][1] .. "," .. b[1][2] .. ") ("
+ .. b[2][1] .. "," .. b[2][2] .. ") ("
+ .. b[1][1] .. "," .. b[1][2] .. ") ("
+ .. b[2][1] .. "," .. b[2][2] .. ") }"
+ else -- notation = tikz
+ bezierstring = "(" .. b[1][1] .. "," .. b[1][2] ..")"
+ for i = 2, #b do
+ bezierstring = bezierstring ..
+ " -- (" .. b[i][1] .. "," .. b[i][2] ..")"
+ end
+ end
+ end
+ else
+ if rev then
+ bezierstring = "(" .. b[#b][#b[#b]-1] .. ","
+ .. b[#b][#b[#b]] ..")" -- initial point
+ for i = #b, 2, -1 do
+ if #b[i] >= 6 then -- cubic bezier spline
+ bezierstring = bezierstring .. " .. controls ("
+ .. b[i][3] .. "," .. b[i][4] ..") and ("
+ .. b[i][1] .. "," .. b[i][2] .. ") .. ("
+ .. b[i-1][#b[i-1]-1] .. "," .. b[i-1][#b[i-1]]..")"
+ else
+ bezierstring = bezierstring .. " ("
+ .. b[i-1][#b[i-1]-1] .. "," .. b[i-1][#b[i-1]] ..")"
+ end
+ end
+ else
+ if notation == "pgfplots" then
+ bezierstring = "\\addplot coordinates {"
+ for i = 1, #b-1 do
+ if #b[i+1] >= 6 then -- cubic bezier spline
+ bezierstring = bezierstring .. "("
+ .. b[i][#b[i]-1] .. "," .. b[i][#b[i]] .. ") ("
+ .. b[i+1][5] .. "," .. b[i+1][6] .. ") ("
+ .. b[i+1][1] .. "," .. b[i+1][2] .. ") ("
+ .. b[i+1][3] .. "," .. b[i+1][4] .. ") "
+ end
+ end
+ bezierstring = bezierstring .. "}"
+ else -- notation = tikz
+ bezierstring = "(" .. b[1][1] .. ","
+ .. b[1][2] ..")" -- initial point
+ for i = 2, #b do
+ if #b[i] >= 6 then -- cubic bezier spline
+ bezierstring = bezierstring .. " .. controls ("
+ .. b[i][1] .. "," .. b[i][2] ..") and ("
+ .. b[i][3] .. "," .. b[i][4] .. ") .. ("
+ .. b[i][5] .. "," .. b[i][6]..")"
+ else
+ bezierstring = bezierstring .. " ("
+ .. b[i][1] .. "," .. b[i][2] ..")"
+ end
+ end
+ end
+ end
+ end
+ end
+ return bezierstring
+end
+
+-- take end points of a graph g of the function f
+-- (from indices starti to endi)
+-- without extrema or inflection points inbetween
+-- and try to approximate it with a cubic bezier curve
+-- (round to rndx and rndy when printing)
+-- if maxerror <= 0, the function will not be recursive anymore
+local function graphtobezierapprox(f,g,starti,endi,maxerror)
+ local px = g[starti][1]
+ local py = g[starti][2]
+ local dp = g[starti][3]
+ local sx = g[endi][1]
+ local sy = g[endi][2]
+ local ds = g[endi][3]
+ -- we compute the corner point c, where the controls would meet
+ local cx = ((dp * px) - (ds * sx) - py + sy) / (dp - ds)
+ local cy = (dp * ((ds * px) - (ds * sx) - py + sy) / (dp - ds)) + py
+ -- now we slide q between p and c & r between s and c
+ -- and search for the best qx and best rx
+ local qx = px+.01*(cx-px)
+ local qy = py+.01*(cy-py)
+ local rx = sx+.01*(cx-sx)
+ local ry = sy+.01*(cy-sy)
+ local err = squareerror(f,g,starti,endi,qx,qy,rx,ry)
+ for i = 2, 99 do
+ for j = 2, 99 do
+ xa = px+i*.01*(cx-px)
+ ya = py+i*.01*(cy-py)
+ xb = sx+j*.01*(cx-sx)
+ yb = sy+j*.01*(cy-sy)
+ -- now check, if xa and xb fit better
+ -- than the last qx and rx did
+ -- (sum of squares must be smaller)
+ local newerror = squareerror(f,g,starti,endi,xa,ya,xb,yb)
+ if newerror < err then
+ qx = xa
+ qy = ya
+ rx = xb
+ ry = yb
+ err = newerror
+ end
+ end
+ end
+ if maxerror > 0 then
+ -- check if it is close enough: (recycling err, xa, ya)
+ err = 0
+ for t = .1, .9, .1 do
+ xa = (1-t)^3*g[starti][1]+3*t*(1-t)^2*qx+3*t^2*(1-t)*rx+t^3*g[endi][1]
+ ya = (1-t)^3*g[starti][2]+3*t*(1-t)^2*qy+3*t^2*(1-t)*ry+t^3*g[endi][2]
+ if abs(ya-f(xa)) > err then
+ err = abs(ya-f(xa))
+ err = abs(ya-f(xa))
+ end
+ end
+ if (err <= maxerror)
+ and qx > -math.huge and qx < math.huge
+ and qy > -math.huge and qy < math.huge
+ and rx > -math.huge and ry < math.huge
+ and sx > -math.huge and sy < math.huge
+ then
+ return {qx,qy,rx,ry,sx,sy}
+ else
+ -- search for an intermediate point where the graph has the same
+ -- slope as the line from the start point to the end point:
+ local interindex = math.floor(.5*starti+.5*endi) -- will change
+ for i = starti + 1, endi - 1 do
+ if abs(g[i][3]-(g[endi][2]-g[starti][2])
+ /(g[endi][1]-g[starti][1]))
+ < abs(g[interindex][3]-(g[endi][2]-g[starti][2])
+ /(g[endi][1]-g[starti][1])) then
+ interindex = i
+ end
+ end
+ local left = graphtobezierapprox(f,g,starti,interindex,maxerror)
+ local right = graphtobezierapprox(f,g,interindex,endi,maxerror)
+ for i=1, #right do --now append the right to the left:
+ left[#left+1] = right[i]
+ end
+ return left
+ end
+ elseif qx > -math.huge and qx < math.huge
+ and qy > -math.huge and qy < math.huge
+ and rx > -math.huge and ry < math.huge
+ and sx > -math.huge and sy < math.huge
+ then
+ return {qx,qy,rx,ry,sx,sy}
+ end
+end
+
+-- like above but exact for quadratic and cubic (if not inverse)
+-- resp. exact for squareroot and cuberoot (if inverse)
+local function graphtobezier(g,starti,endi,isinverse)
+ local px = g[starti][1]
+ local py = g[starti][2]
+ local dp = g[starti][3]
+ local sx = g[endi][1]
+ local sy = g[endi][2]
+ local ds = g[endi][3]
+ local qx = px+(sx-px)/3
+ local rx = px+2*(sx-px)/3
+ local qy = py+(qx-px)*dp
+ local ry = sy+(rx-sx)*ds
+ if isinverse then
+ return {qy,qx,ry,rx,sy,sx}
+ else
+ return {qx,qy,rx,ry,sx,sy}
+ end
+end
+
+-- main function
+function bezierplot(functionstring,xminstring,xmaxstring,yminstring,ymaxstring,samplesstring,notation)
+ local fstringreplaced = string.gsub(functionstring, "%*%*", "^")
+ local f = assert(load("local x = ...; return " .. fstringreplaced))
+ local xmin = evaluate(xminstring)
+ local xmax = evaluate(xmaxstring)
+ local ymin = evaluate(yminstring)
+ local ymax = evaluate(ymaxstring)
+ local samples = evaluate(samplesstring)
+ local isreverse = false
+ if xmin > xmax then
+ isreverse = true
+ elseif xmin == xmax then
+ xmax = xmin + 10
+ end
+ xmin, xmax = math.min(xmin,xmax), math.max(xmin,xmax)
+ if ymin == ymax then
+ ymax = ymin + 10
+ end
+ ymin, ymax = math.min(ymin,ymax), math.max(ymin,ymax)
+ local xsteps = 50000
+ -- if samples < 2 the samples will be chosen as wisely as possible
+ local arbitrary_samples = true
+ if samples >= 2 then
+ arbitrary_samples = false
+ xsteps = (samples-1)*math.max(2,math.floor(xsteps/(samples-1)))
+ end
+ local xstep = (xmax-xmin)/xsteps
+ -- the output of the x coordinates will be rounded to rndx digits
+ local rndx = math.max(0,math.floor(5.5-log(xmax-xmin)/log(10)))
+ local xerror = abs(xmax-xmin)/(10^rndx)
+ -- the output of the y coordinates will be rounded to rndy digits
+ local rndy = math.max(0,math.floor(5.5-log(ymax-ymin)/log(10)))
+ local yerror = (ymax-ymin)/(10^rndy)
+ -- determine parts of the graph that are inside window
+ local graphs = {} -- graph split to the connected parts
+ local graph = {} -- graphs concatenated (needed for function type)
+ local outside = true -- value is outside window
+ local i = 0
+ local j = 0
+ local yminreal -- determine the real minimimum of the y coord.
+ local ymaxreal -- just decring
+ local yminrealfound = false
+ local ymaxrealfound = false
+ for n = 0, xsteps do
+ local x = xmin + n/xsteps*(xmax-xmin)
+ if n == xsteps then
+ x = xmax
+ end
+ local y = f(x)
+ if (y >= ymin-.1*yerror and ymin ~= -huge or y > ymin and ymin == -huge)
+ and (y <= ymax+.1*yerror and ymax ~= huge or y < ymax and ymax == huge)
+ then -- inside
+ if outside then -- if it was outside before
+ outside = false
+ j = 0
+ i = i + 1
+ graphs[i] = {}
+ end
+ j = j + 1
+ graphs[i][j] = {x,y}
+ graph[#graph+1] = {x,y}
+ if not yminrealfound or yminrealfound and y < yminreal then
+ yminreal = y
+ yminrealfound = true
+ end
+ if not ymaxrealfound or ymaxrealfound and y > ymaxreal then
+ ymaxreal = y
+ ymaxrealfound = true
+ end
+ else
+ outside = true
+ end
+ end
+
+ -- some redefinitions
+ if #graph ~= 0 and yminreal ~= ymaxreal then
+ ymin = yminreal
+ ymax = ymaxreal
+ rndy = math.max(0,math.floor(5.5-log(ymax-ymin)/log(10)))
+ yerror = (ymax-ymin)/(10^rndy)
+ end
+
+ -- check for the function type (for this, we need the concatenated
+ -- parts of the graph)
+ -- go through the connected parts
+ local functiontype = "unknown"
+ local a, b, c, d -- possible function parameter
+ -- check for affine functions:
+ local l = #graph
+ a, b = parameters_affine(graph[1][1],graph[1][2],
+ graph[l][1],graph[l][2])
+ if do_parameters_fit(a,b,0,0,"a*x+b",graph,yerror,false) then
+ functiontype = "affine"
+ else -- check for cubic functions (includes quadratic functions)
+ a, b, c, d = parameters_cubic(graph[1][1],graph[1][2],
+ graph[math.floor(l/3)][1],graph[math.floor(l/3)][2],
+ graph[math.floor(2*l/3)][1],graph[math.floor(2*l/3)][2],
+ graph[l][1],graph[l][2])
+ if do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d",graph,
+ yerror,false) then
+ functiontype = "cubic"
+ else -- check for cuberoot functions (includes squareroots)
+ a, b, c, d = parameters_cubic(graph[1][2],graph[1][1],
+ graph[math.floor(l/3)][2],graph[math.floor(l/3)][1],
+ graph[math.floor(2*l/3)][2],graph[math.floor(2*l/3)][1],
+ graph[l][2],graph[l][1])
+ if do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d",graph,
+ xerror,true) then
+ functiontype = "cuberoot"
+ end
+ end
+ end
+
+ local bezierpoints = {}
+ -- the bezier path (0,1) .. controls
+ -- (2,3) and (4,5) .. (6,7) .. controls
+ -- (8,9) and (10,11) .. (12,13)
+ -- will be stored as
+ -- bezierpoints={{0,1},{2,3,4,5,6,7},{8,9,10,11,12,13}}
+
+ if functiontype == "affine" then
+ if arbitrary_samples then
+ bezierpoints = {{graph[1][1],graph[1][2]},{graph[#graph][1],
+ graph[#graph][2]}}
+ else -- we can here savely assume that graphs has only one part,
+ -- therefore graphs[1]=graph
+ for i = 1, #graph, math.floor(xsteps/(samples-1)) do
+ bezierpoints[#bezierpoints+1] = {graph[i][1],graph[i][2]}
+ end
+ end
+ elseif functiontype == "cubic" then
+ local extrema_inflections = {} -- store the extrema and
+ -- inflection points for arbitrary samples
+ if arbitrary_samples then
+ if math.abs(a) < yerror*1e-10 then -- quadratic case (one extremum)
+ if is_in_window(-c/(2*b),(-c^2+4*b*d)/(4*b),xmin,xmax,
+ ymin,ymax) then
+ extrema_inflections = {{-c/(2*b),(-c^2+4*b*d)/(4*b)}}
+ end
+ else -- cubic case (two extrema and one inflection point)
+ -- we order the points with the help of sgn
+ -- check for first extrema
+ if is_in_window((-sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a),
+ (2*b^3+27*a^2*d-9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)*
+ (2*b^2-6*a*c))/(27*a^2),xmin,xmax,ymin,ymax) then
+ extrema_inflections[#extrema_inflections+1] =
+ {(-sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a),(2*b^3+27*a^2*d-
+ 9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)*(2*b^2-6*a*c))/(27*a^2)}
+ end
+ -- check for inflection point (has to be inbetween)
+ if is_in_window(-b/(3*a),(2*b^3+27*a^2*d-9*a*b*c)
+ /(27*a^2),xmin,xmax,ymin,ymax) then
+ extrema_inflections[#extrema_inflections+1]={-b/(3*a),
+ (2*b^3+27*a^2*d-9*a*b*c)/(27*a^2)}
+ end
+ -- check for second extrema
+ if is_in_window((sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a),
+ (2*b^3+27*a^2*d-9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)*
+ (-2*b^2+6*a*c))/(27*a^2),xmin,xmax,ymin,ymax) then
+ extrema_inflections[#extrema_inflections+1] =
+ {(sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a),(2*b^3+27*a^2*d-
+ 9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)*(-2*b^2+6*a*c))/(27*a^2)}
+ end
+ end
+ end
+ for part = 1, #graphs do
+ bezierpoints[#bezierpoints+1] = {graphs[part][1][1],
+ graphs[part][1][2]} -- initial points
+ local graphsamples = {}-- will be the graph reduced to the
+ -- samples (or the most important points)
+ local dg -- will be the differentiated graph
+ if arbitrary_samples then -- add extrema and inflection
+ -- points to the border points
+ graphsamples = {{graphs[part][1][1],
+ graphs[part][1][2]}}
+ for j = 1, #extrema_inflections do
+ if extrema_inflections[j][1] > math.min(
+ graphs[part][1][1] ,graphs[part][#graphs[part]][1])
+ and extrema_inflections[j][1] < math.max(
+ graphs[part][1][1] ,graphs[part][#graphs[part]][1])
+ then
+ graphsamples[#graphsamples+1] =
+ {extrema_inflections[j][1],
+ extrema_inflections[j][2]}
+ end
+ end
+ graphsamples[#graphsamples+1] =
+ {graphs[part][#graphs[part]][1],
+ graphs[part][#graphs[part]][2]}
+ else
+ for i = 1, #graphs[part], xsteps/(samples-1) do
+ graphsamples[#graphsamples+1] =
+ {graphs[part][i][1],graphs[part][i][2]}
+ end
+ end
+ dg = diffgraphcubic(graphsamples,a,b,c,d,false)
+ for i = 2, #dg do
+ bezierpoints[#bezierpoints+1] = graphtobezier(dg,i-1,i,false)
+ end
+ end
+ elseif functiontype == "cuberoot" then
+ local inflection = {} -- store the inflection point
+ if arbitrary_samples and math.abs(a) ~= 0
+ and is_in_window((2*b^3+27*a^2*d-9*a*b*c)/(27*a^2),-b/(3*a),
+ xmin,xmax,ymin,ymax) then
+ inflection = {(2*b^3+27*a^2*d-9*a*b*c)/(27*a^2),-b/(3*a)}
+ end
+ -- (there cannot be more than one part)
+ bezierpoints[#bezierpoints+1] = {graphs[1][1][1],
+ graphs[1][1][2]} -- initial points
+ local graphsamples = {}-- will be the graph reduced to the
+ -- samples (or the most important points)
+ local dg -- will be the differentiated graph
+ if arbitrary_samples then -- add inflection point (if exis.)
+ graphsamples = {{graphs[1][1][1],
+ graphs[1][1][2]}}
+ if #inflection > 0 and inflection[1] > math.min(
+ graphs[1][1][1],graphs[1][#graphs[1]][1])
+ and inflection[1] < math.max(
+ graphs[1][1][1],graphs[1][#graphs[1]][1])
+ then
+ graphsamples[#graphsamples+1] =
+ {inflection[1],inflection[2]}
+ end
+ graphsamples[#graphsamples+1] =
+ {graphs[1][#graphs[1]][1],
+ graphs[1][#graphs[1]][2]}
+ else
+ for i = 1, #graphs[1], xsteps/(samples-1) do
+ graphsamples[#graphsamples+1] =
+ {graphs[1][i][1],graphs[1][i][2]}
+ end
+ end
+ dg = diffgraphcubic(graphsamples,a,b,c,d,true)
+ for i = 2, #dg do
+ bezierpoints[#bezierpoints+1] = graphtobezier(dg,i-1,i,true)
+ end
+ else
+ ---------- generic case (no special function) ----------------
+ if arbitrary_samples then
+ -- go through the connected parts
+ for part = 1, #graphs do
+ local dg = diffgraph(f,graphs[part],xstep)
+ --printdifftable(dg) -- for debugging
+ bezierpoints[#bezierpoints+1] = {dg[1][1],dg[1][2]}
+ local startindex = 1
+ for k = 2, #dg do
+ if dg[k][5] or dg[k][6] then -- extrema and inflection points
+ local tobeadded = graphtobezierapprox(
+ f,dg,startindex,k,10*yerror)
+ -- tobeadded may contain a multiple of 6 entries
+ -- e.g. {1,2,3,4,5,6,7,8,9,10,11,12}
+ for i = 1, math.floor(#tobeadded/6) do
+ bezierpoints[#bezierpoints+1] = {}
+ for j = 1, 6 do
+ bezierpoints[#bezierpoints][j] = tobeadded[(i-1)*6+j]
+ end
+ end
+ startindex = k
+ end
+ end
+ if startindex ~= #dg then -- if no special points inbetween
+ local tobeadded = graphtobezierapprox(f,dg,
+ startindex,#dg,10*yerror)
+ -- tobeadded may contain a multiple of 6 entries
+ -- e.g. {1,2,3,4,5,6,7,8,9,10,11,12}
+ for i = 1, math.floor(#tobeadded/6) do
+ bezierpoints[#bezierpoints+1] = {}
+ for j = 1, 6 do
+ bezierpoints[#bezierpoints][j] = tobeadded[(i-1)*6+j]
+ end
+ end
+ end
+ end
+ else -- fixed samples in the generic case
+ -- go through the connected parts
+ for part = 1, #graphs do
+ local dg = diffgraphsimple(f,graphs[part],xstep,
+ math.floor(0.5+xsteps/(samples-1)))
+ bezierpoints[#bezierpoints+1] = {dg[1][1],dg[1][2]} -- initial points
+ for i = 2, #dg do
+ bezierpoints[#bezierpoints+1] = graphtobezier(dg,i-1,i,false)
+ end
+ end
+ end
+ end
+ return beziertabletostring(bezierpoints,rndx,rndy,isreverse,notation)
+end
+
+-- main program --
+
+if not pcall(debug.getlocal, 4, 1) then
+--if debug.getinfo(3) == nil then
+ if #arg >= 1 then
+ local xmin = -5
+ local xmax = 5
+ if #arg >= 2 then
+ xmin = arg[2]
+ end
+ if #arg >= 3 then
+ xmax = arg[3]
+ end
+ local ymin = -5
+ local ymax = 5
+ if #arg >= 4 then
+ ymin = arg[4]
+ end
+ if #arg >= 5 then
+ ymax = arg[5]
+ end
+ local samples = 0
+ if #arg >= 6 then
+ samples = arg[6]
+ end
+ local notation = "tikz"
+ if #arg >= 7 then
+ notation = arg[7]
+ end
+ print(bezierplot(arg[1],xmin,xmax,ymin,ymax,samples,notation))
+ end
+end
+
+
+
diff --git a/macros/luatex/latex/bezierplot/bezierplot.sty b/macros/luatex/latex/bezierplot/bezierplot.sty
new file mode 100644
index 0000000000..94be30b4dc
--- /dev/null
+++ b/macros/luatex/latex/bezierplot/bezierplot.sty
@@ -0,0 +1,17 @@
+\NeedsTeXFormat{LaTeX2e}
+\ProvidesPackage{bezierplot}[2019/06/18 bezierplot]
+\RequirePackage{xparse}
+\RequirePackage{iftex}
+\ifLuaTeX
+ \directlua{require("bezierplot")}
+ \DeclareExpandableDocumentCommand{\xbezierplot}{O{-5} O{5} O{-5} O{5} O{0} O{tikz} m}{%
+ \directlua{tex.sprint(bezierplot("#7","#1","#2","#3","#4","#5","#6"))}
+ }
+\else
+ \let\xpandblinpt\@@input
+ \DeclareExpandableDocumentCommand{\xbezierplot}{O{-5} O{5} O{-5} O{5} O{0} O{tikz} m}{%
+ \xpandblinpt|"bezierplot '#7' '#1' '#2' '#3' '#4' '#5' '#6'"
+ }
+\fi
+\providecommand\bezierplot{\romannumeral`\^^@\xbezierplot}
+\endinput \ No newline at end of file