Best piecewise with $n$ points for a portion of a log function. Let $f(x) = |\log_2(x)|$ for $x$ belonging to the domain $(0,1]$. I would like to know if there is an algorithm to fit $f(x)$ using a piecewise linear function g(x) in an optimal way?. That is, on input $n$, the algorithm needs to output the better linear approximation using $n$ pieces. Also $g(x)<=f(x)$. Maybe there is a function/procedure in Python or other programming language?

You can always consider the norm
$$\Phi(\alpha ,\beta)=\int_a^b \Bigg[\frac{\log (x)}{\log (2)}-(\alpha +\beta  x)\Bigg]^2$$ and minimize it with respect to $(\alpha ,\beta)$. This is equivalent to the sum of squares of a linear regression for an infinite number of points for $x \in
(\alpha ,\beta)$.
The antiderivative is quite simple and writing
$$\frac{\partial \Phi(\alpha ,\beta)}{\partial \alpha}=\frac{\partial \Phi(\alpha ,\beta)}{\partial \beta}=0$$ gives two linear equations in $(\alpha ,\beta)$.
Skipping the intermediate calculations, you should get
$$ \beta=\frac{3}{(a-b)^3\log (2) }\Bigg[(a^2-b^2) +2 a b \log \left(\frac{b}{a}\right)\Bigg]$$
$$\alpha=\frac{(b-a) (\beta  \log (2) (a+b)+2)+2 a \log (a)-2 b \log (b)}{2  (a-b) \log (2) }$$
Try for $a=\frac 12$ and $b=\frac 12+\frac 1{10}=\frac 35$ which give
$$\alpha=\frac{-365-1992 \log \left(\frac{5}{3}\right)+1990 \log (2)}{2 \log (2)}\qquad\qquad
\beta=\frac{330-3600 \coth ^{-1}(11)}{\log (2)}$$ gives a maximum absolute error of $0.004$ at the bounds.
But, making $a=\frac 12$ and $b=\frac 12+\frac 1{100}=\frac {51}{100}$, the maximum absolute error is $0.00005$ at the bounds.
But, because of the curvatur of the logarithmic function, I do not see how your constraint $g(x) \leq f(x)$ could (at least easily) could be satisfied.
Edit
Making $b=a+\epsilon$ and using series expansions with $t=\frac \epsilon a$
$$\alpha=\frac 1{2\log(2)}\Bigg[2 \log \left(\frac{a}{e}\right)+\sum_{n=1}^\infty (-1)^{n-1}\, c_n\,t^n\Bigg]$$ where the $c_n$ form the sequence
$$\left\{1,\frac{13}{30},\frac{4}{15},\frac{13}{70},\frac{29}{210},\frac{3}{28},\frac{3}{35},\frac{139}{1980},\cdots\right\}$$
$$\beta=\frac{1}{a \log (2)}\sum_{n=0}^\infty (-1)^{n}\, d_n\,t^n\Bigg]$$where the $d_n$ form the sequence
$$\left\{1,\frac{1}{2},\frac{3}{10},\frac{1}{5},\frac{1}{7},\frac{3}{28},\frac{1}{12},\frac{1}{15},\frac{3}{55},\cdots\right\}$$ and
$$\Phi_{\text{min}}\sim\frac{\epsilon ^5}{720\, a^4 \log ^2(2)}$$
