Skip to content
Snippets Groups Projects
Verified Commit 1c1fbb05 authored by Herbers, Maik's avatar Herbers, Maik
Browse files

import

parents
No related branches found
No related tags found
No related merge requests found
(authorizations
(version 0)
(("D5BA 4708 FA7D 0AFD 9C0E 6556 ECD7 F82F 5327 404C" ; signing subkey
(name "Maik Herbers"))))
guix.scm 0 → 100644
;;; SPDX-License-Identifier: GPL-3.0-or-later
(use-modules (gnu packages image)
(gnu packages multiprecision)
(gnu packages pkg-config)
(guix build-system meson)
(guix gexp)
(guix licenses)
(guix packages)
(ice-9 popen)
(ice-9 rdelim))
(define polyplot
(let*
((%source-dir (dirname (current-filename)))
(cmd (string-append "git -C " %source-dir " describe --always"))
(port (open-pipe cmd OPEN_READ))
(%git-version (read-line port))
(select? (lambda (name _)
(not (string=? name
(string-append %source-dir "/.git"))))))
(close port)
(package
(name "polyplot")
(version %git-version)
(source (local-file %source-dir
#:recursive? #t
#:select? select?))
(build-system meson-build-system)
(arguments (list #:tests? #f))
(inputs (list libpng mpc mpfr))
(native-inputs (list pkg-config))
(synopsis "Plotter for complex power series")
(description "This package creates plots of power series on their domain of
convergence. It uses the perceptually uniform colormap @code{cividis}
to represent phase.")
(home-page (string-append "file://localhost" %source-dir))
(license gpl3+))))
polyplot
/* SPDX-License-Identifier: GPL-3.0-or-later */
#pragma once
#include "common.h"
/*
* Generated from matplotlib by
* import matplotlib as mpl
* cmap = mpl.cm.get_cmap('cividis')
* to_hex = lambda val: f"0x{int(round(255 * val)):02x}"
* for i in range(cmap.N):
* print("{{ {}, {}, {}, {} }},".format(*map(to_hex, cmap(i))))
*/
struct rgba_pix cividis_lut[] = {
{ 0x00, 0x22, 0x4e, 0xff },
{ 0x00, 0x23, 0x4f, 0xff },
{ 0x00, 0x24, 0x51, 0xff },
{ 0x00, 0x25, 0x53, 0xff },
{ 0x00, 0x25, 0x54, 0xff },
{ 0x00, 0x26, 0x56, 0xff },
{ 0x00, 0x27, 0x58, 0xff },
{ 0x00, 0x28, 0x59, 0xff },
{ 0x00, 0x28, 0x5b, 0xff },
{ 0x00, 0x29, 0x5d, 0xff },
{ 0x00, 0x2a, 0x5f, 0xff },
{ 0x00, 0x2a, 0x61, 0xff },
{ 0x00, 0x2b, 0x62, 0xff },
{ 0x00, 0x2c, 0x64, 0xff },
{ 0x00, 0x2c, 0x66, 0xff },
{ 0x00, 0x2d, 0x68, 0xff },
{ 0x00, 0x2e, 0x6a, 0xff },
{ 0x00, 0x2e, 0x6c, 0xff },
{ 0x00, 0x2f, 0x6d, 0xff },
{ 0x00, 0x30, 0x6f, 0xff },
{ 0x00, 0x30, 0x70, 0xff },
{ 0x00, 0x31, 0x70, 0xff },
{ 0x00, 0x31, 0x71, 0xff },
{ 0x01, 0x32, 0x71, 0xff },
{ 0x05, 0x33, 0x71, 0xff },
{ 0x08, 0x33, 0x70, 0xff },
{ 0x0c, 0x34, 0x70, 0xff },
{ 0x0f, 0x35, 0x70, 0xff },
{ 0x12, 0x35, 0x70, 0xff },
{ 0x14, 0x36, 0x70, 0xff },
{ 0x16, 0x37, 0x70, 0xff },
{ 0x18, 0x37, 0x6f, 0xff },
{ 0x1a, 0x38, 0x6f, 0xff },
{ 0x1c, 0x39, 0x6f, 0xff },
{ 0x1e, 0x3a, 0x6f, 0xff },
{ 0x20, 0x3a, 0x6f, 0xff },
{ 0x21, 0x3b, 0x6e, 0xff },
{ 0x23, 0x3c, 0x6e, 0xff },
{ 0x24, 0x3c, 0x6e, 0xff },
{ 0x26, 0x3d, 0x6e, 0xff },
{ 0x27, 0x3e, 0x6e, 0xff },
{ 0x29, 0x3f, 0x6e, 0xff },
{ 0x2a, 0x3f, 0x6d, 0xff },
{ 0x2b, 0x40, 0x6d, 0xff },
{ 0x2d, 0x41, 0x6d, 0xff },
{ 0x2e, 0x41, 0x6d, 0xff },
{ 0x2f, 0x42, 0x6d, 0xff },
{ 0x31, 0x43, 0x6d, 0xff },
{ 0x32, 0x43, 0x6d, 0xff },
{ 0x33, 0x44, 0x6d, 0xff },
{ 0x34, 0x45, 0x6c, 0xff },
{ 0x35, 0x45, 0x6c, 0xff },
{ 0x36, 0x46, 0x6c, 0xff },
{ 0x38, 0x47, 0x6c, 0xff },
{ 0x39, 0x48, 0x6c, 0xff },
{ 0x3a, 0x48, 0x6c, 0xff },
{ 0x3b, 0x49, 0x6c, 0xff },
{ 0x3c, 0x4a, 0x6c, 0xff },
{ 0x3d, 0x4a, 0x6c, 0xff },
{ 0x3e, 0x4b, 0x6c, 0xff },
{ 0x3f, 0x4c, 0x6c, 0xff },
{ 0x40, 0x4c, 0x6c, 0xff },
{ 0x41, 0x4d, 0x6c, 0xff },
{ 0x42, 0x4e, 0x6c, 0xff },
{ 0x43, 0x4e, 0x6c, 0xff },
{ 0x44, 0x4f, 0x6c, 0xff },
{ 0x45, 0x50, 0x6c, 0xff },
{ 0x46, 0x51, 0x6c, 0xff },
{ 0x47, 0x51, 0x6c, 0xff },
{ 0x48, 0x52, 0x6c, 0xff },
{ 0x49, 0x53, 0x6c, 0xff },
{ 0x4a, 0x53, 0x6c, 0xff },
{ 0x4b, 0x54, 0x6c, 0xff },
{ 0x4c, 0x55, 0x6c, 0xff },
{ 0x4d, 0x55, 0x6c, 0xff },
{ 0x4e, 0x56, 0x6c, 0xff },
{ 0x4f, 0x57, 0x6c, 0xff },
{ 0x50, 0x57, 0x6c, 0xff },
{ 0x51, 0x58, 0x6d, 0xff },
{ 0x52, 0x59, 0x6d, 0xff },
{ 0x53, 0x5a, 0x6d, 0xff },
{ 0x54, 0x5a, 0x6d, 0xff },
{ 0x55, 0x5b, 0x6d, 0xff },
{ 0x55, 0x5c, 0x6d, 0xff },
{ 0x56, 0x5c, 0x6d, 0xff },
{ 0x57, 0x5d, 0x6d, 0xff },
{ 0x58, 0x5e, 0x6d, 0xff },
{ 0x59, 0x5e, 0x6e, 0xff },
{ 0x5a, 0x5f, 0x6e, 0xff },
{ 0x5b, 0x60, 0x6e, 0xff },
{ 0x5c, 0x61, 0x6e, 0xff },
{ 0x5d, 0x61, 0x6e, 0xff },
{ 0x5e, 0x62, 0x6e, 0xff },
{ 0x5e, 0x63, 0x6f, 0xff },
{ 0x5f, 0x63, 0x6f, 0xff },
{ 0x60, 0x64, 0x6f, 0xff },
{ 0x61, 0x65, 0x6f, 0xff },
{ 0x62, 0x65, 0x6f, 0xff },
{ 0x63, 0x66, 0x70, 0xff },
{ 0x64, 0x67, 0x70, 0xff },
{ 0x65, 0x68, 0x70, 0xff },
{ 0x65, 0x68, 0x70, 0xff },
{ 0x66, 0x69, 0x70, 0xff },
{ 0x67, 0x6a, 0x71, 0xff },
{ 0x68, 0x6a, 0x71, 0xff },
{ 0x69, 0x6b, 0x71, 0xff },
{ 0x6a, 0x6c, 0x71, 0xff },
{ 0x6b, 0x6d, 0x72, 0xff },
{ 0x6c, 0x6d, 0x72, 0xff },
{ 0x6c, 0x6e, 0x72, 0xff },
{ 0x6d, 0x6f, 0x72, 0xff },
{ 0x6e, 0x6f, 0x73, 0xff },
{ 0x6f, 0x70, 0x73, 0xff },
{ 0x70, 0x71, 0x73, 0xff },
{ 0x71, 0x72, 0x74, 0xff },
{ 0x72, 0x72, 0x74, 0xff },
{ 0x72, 0x73, 0x74, 0xff },
{ 0x73, 0x74, 0x75, 0xff },
{ 0x74, 0x74, 0x75, 0xff },
{ 0x75, 0x75, 0x75, 0xff },
{ 0x76, 0x76, 0x76, 0xff },
{ 0x77, 0x77, 0x76, 0xff },
{ 0x77, 0x77, 0x77, 0xff },
{ 0x78, 0x78, 0x77, 0xff },
{ 0x79, 0x79, 0x77, 0xff },
{ 0x7a, 0x7a, 0x78, 0xff },
{ 0x7b, 0x7a, 0x78, 0xff },
{ 0x7c, 0x7b, 0x78, 0xff },
{ 0x7d, 0x7c, 0x78, 0xff },
{ 0x7e, 0x7c, 0x78, 0xff },
{ 0x7e, 0x7d, 0x78, 0xff },
{ 0x7f, 0x7e, 0x78, 0xff },
{ 0x80, 0x7f, 0x78, 0xff },
{ 0x81, 0x7f, 0x78, 0xff },
{ 0x82, 0x80, 0x79, 0xff },
{ 0x83, 0x81, 0x79, 0xff },
{ 0x84, 0x82, 0x79, 0xff },
{ 0x85, 0x82, 0x79, 0xff },
{ 0x86, 0x83, 0x79, 0xff },
{ 0x87, 0x84, 0x78, 0xff },
{ 0x88, 0x85, 0x78, 0xff },
{ 0x89, 0x85, 0x78, 0xff },
{ 0x8a, 0x86, 0x78, 0xff },
{ 0x8b, 0x87, 0x78, 0xff },
{ 0x8c, 0x88, 0x78, 0xff },
{ 0x8d, 0x88, 0x78, 0xff },
{ 0x8e, 0x89, 0x78, 0xff },
{ 0x8f, 0x8a, 0x78, 0xff },
{ 0x90, 0x8b, 0x78, 0xff },
{ 0x91, 0x8b, 0x78, 0xff },
{ 0x92, 0x8c, 0x78, 0xff },
{ 0x92, 0x8d, 0x78, 0xff },
{ 0x93, 0x8e, 0x78, 0xff },
{ 0x94, 0x8e, 0x77, 0xff },
{ 0x95, 0x8f, 0x77, 0xff },
{ 0x96, 0x90, 0x77, 0xff },
{ 0x97, 0x91, 0x77, 0xff },
{ 0x98, 0x92, 0x77, 0xff },
{ 0x99, 0x92, 0x77, 0xff },
{ 0x9a, 0x93, 0x76, 0xff },
{ 0x9b, 0x94, 0x76, 0xff },
{ 0x9c, 0x95, 0x76, 0xff },
{ 0x9d, 0x95, 0x76, 0xff },
{ 0x9e, 0x96, 0x76, 0xff },
{ 0x9f, 0x97, 0x75, 0xff },
{ 0xa0, 0x98, 0x75, 0xff },
{ 0xa1, 0x99, 0x75, 0xff },
{ 0xa2, 0x99, 0x75, 0xff },
{ 0xa3, 0x9a, 0x74, 0xff },
{ 0xa4, 0x9b, 0x74, 0xff },
{ 0xa5, 0x9c, 0x74, 0xff },
{ 0xa6, 0x9c, 0x74, 0xff },
{ 0xa7, 0x9d, 0x73, 0xff },
{ 0xa8, 0x9e, 0x73, 0xff },
{ 0xa9, 0x9f, 0x73, 0xff },
{ 0xaa, 0xa0, 0x73, 0xff },
{ 0xab, 0xa0, 0x72, 0xff },
{ 0xac, 0xa1, 0x72, 0xff },
{ 0xad, 0xa2, 0x72, 0xff },
{ 0xae, 0xa3, 0x71, 0xff },
{ 0xaf, 0xa4, 0x71, 0xff },
{ 0xb0, 0xa5, 0x71, 0xff },
{ 0xb1, 0xa5, 0x70, 0xff },
{ 0xb3, 0xa6, 0x70, 0xff },
{ 0xb4, 0xa7, 0x6f, 0xff },
{ 0xb5, 0xa8, 0x6f, 0xff },
{ 0xb6, 0xa9, 0x6f, 0xff },
{ 0xb7, 0xa9, 0x6e, 0xff },
{ 0xb8, 0xaa, 0x6e, 0xff },
{ 0xb9, 0xab, 0x6d, 0xff },
{ 0xba, 0xac, 0x6d, 0xff },
{ 0xbb, 0xad, 0x6d, 0xff },
{ 0xbc, 0xae, 0x6c, 0xff },
{ 0xbd, 0xae, 0x6c, 0xff },
{ 0xbe, 0xaf, 0x6b, 0xff },
{ 0xbf, 0xb0, 0x6b, 0xff },
{ 0xc0, 0xb1, 0x6a, 0xff },
{ 0xc1, 0xb2, 0x6a, 0xff },
{ 0xc2, 0xb3, 0x69, 0xff },
{ 0xc3, 0xb3, 0x69, 0xff },
{ 0xc4, 0xb4, 0x68, 0xff },
{ 0xc5, 0xb5, 0x68, 0xff },
{ 0xc6, 0xb6, 0x67, 0xff },
{ 0xc7, 0xb7, 0x67, 0xff },
{ 0xc8, 0xb8, 0x66, 0xff },
{ 0xc9, 0xb9, 0x65, 0xff },
{ 0xcb, 0xb9, 0x65, 0xff },
{ 0xcc, 0xba, 0x64, 0xff },
{ 0xcd, 0xbb, 0x63, 0xff },
{ 0xce, 0xbc, 0x63, 0xff },
{ 0xcf, 0xbd, 0x62, 0xff },
{ 0xd0, 0xbe, 0x62, 0xff },
{ 0xd1, 0xbf, 0x61, 0xff },
{ 0xd2, 0xc0, 0x60, 0xff },
{ 0xd3, 0xc0, 0x5f, 0xff },
{ 0xd4, 0xc1, 0x5f, 0xff },
{ 0xd5, 0xc2, 0x5e, 0xff },
{ 0xd6, 0xc3, 0x5d, 0xff },
{ 0xd7, 0xc4, 0x5c, 0xff },
{ 0xd9, 0xc5, 0x5c, 0xff },
{ 0xda, 0xc6, 0x5b, 0xff },
{ 0xdb, 0xc7, 0x5a, 0xff },
{ 0xdc, 0xc8, 0x59, 0xff },
{ 0xdd, 0xc8, 0x58, 0xff },
{ 0xde, 0xc9, 0x58, 0xff },
{ 0xdf, 0xca, 0x57, 0xff },
{ 0xe0, 0xcb, 0x56, 0xff },
{ 0xe1, 0xcc, 0x55, 0xff },
{ 0xe2, 0xcd, 0x54, 0xff },
{ 0xe4, 0xce, 0x53, 0xff },
{ 0xe5, 0xcf, 0x52, 0xff },
{ 0xe6, 0xd0, 0x51, 0xff },
{ 0xe7, 0xd1, 0x50, 0xff },
{ 0xe8, 0xd2, 0x4f, 0xff },
{ 0xe9, 0xd3, 0x4e, 0xff },
{ 0xea, 0xd3, 0x4c, 0xff },
{ 0xeb, 0xd4, 0x4b, 0xff },
{ 0xed, 0xd5, 0x4a, 0xff },
{ 0xee, 0xd6, 0x49, 0xff },
{ 0xef, 0xd7, 0x48, 0xff },
{ 0xf0, 0xd8, 0x46, 0xff },
{ 0xf1, 0xd9, 0x45, 0xff },
{ 0xf2, 0xda, 0x44, 0xff },
{ 0xf3, 0xdb, 0x42, 0xff },
{ 0xf5, 0xdc, 0x41, 0xff },
{ 0xf6, 0xdd, 0x3f, 0xff },
{ 0xf7, 0xde, 0x3e, 0xff },
{ 0xf8, 0xdf, 0x3c, 0xff },
{ 0xf9, 0xe0, 0x3a, 0xff },
{ 0xfb, 0xe1, 0x38, 0xff },
{ 0xfc, 0xe2, 0x36, 0xff },
{ 0xfd, 0xe3, 0x34, 0xff },
{ 0xfe, 0xe4, 0x34, 0xff },
{ 0xfe, 0xe5, 0x35, 0xff },
{ 0xfe, 0xe6, 0x36, 0xff },
{ 0xfe, 0xe8, 0x38, 0xff },
};
/* SPDX-License-Identifier: GPL-3.0-or-later */
#pragma once
#include "common.h"
void get_colour_mono(struct rgba_pix *colour, const mpc_t val, struct temporaries tmps);
void get_colour_cividis(struct rgba_pix *colour, const mpc_t val, struct temporaries tmps);
/* SPDX-License-Identifier: GPL-3.0-or-later */
#pragma once
#include <stdint.h>
#include <mpc.h>
#include <mpfr.h>
struct rgba_pix {
uint8_t r;
uint8_t g;
uint8_t b;
uint8_t a;
};
struct temporaries {
mpfr_t pi;
mpfr_t mpfr;
mpc_t mpc;
};
struct temporaries init_temporaries(uint64_t prec);
void destroy_temporaries(struct temporaries tmps);
/* SPDX-License-Identifier: GPL-3.0-or-later */
#pragma once
#include "common.h"
struct polynomial {
uint64_t degree;
mpc_t *coeffs;
};
void horner_eval(mpc_t res, const mpc_t z, const struct polynomial p);
# SPDX-License-Identifier: GPL-3.0-or-later
configure_file(output : 'config.h', configuration : conf_data)
/* SPDX-License-Identifier: GPL-3.0-or-later */
#pragma once
#include <png.h>
#include "common.h"
struct pic {
png_image im;
struct rgba_pix *pix;
};
struct pic init_pic(uint64_t height, uint64_t width);
void save_to_png(struct pic pic, const char *out);
void free_pic(struct pic pic);
# SPDX-License-Identifier: GPL-3.0-or-later
project(
'polyplot',
'c',
default_options : ['c_std=c18', 'warning_level=3'],
license: 'GPL-3.0-or-later',
version : 'v0.1'
)
add_project_arguments(
[
'-D_POSIX_C_SOURCE=200809L',
'-Wpedantic',
],
language : 'c'
)
mpc = declare_dependency(link_args : '-lmpc')
deps = [
dependency('libpng'),
mpc,
dependency('mpfr'),
dependency('openmp'),
]
conf_data = configuration_data()
conf_data.set_quoted('POLYPLOT_VERSION', meson.project_version())
inc = include_directories('include')
subdir('include')
subdir('src')
/* SPDX-License-Identifier: GPL-3.0-or-later */
#include "colour.h"
#include "cividis_lut.h"
void get_colour_mono(struct rgba_pix *colour, const mpc_t val, struct temporaries tmps)
{
/*
* Use
* 1/π * arctan(log(|val|^α + 1)) ∈ [0, 1]
* as brightness.
*
* See [0] section 2.2.2 for more.
*
* [0]: https://arxiv.org/abs/2002.05234
*/
mpc_abs(tmps.mpfr, val, MPFR_RNDN);
mpfr_rootn_ui(tmps.mpfr, tmps.mpfr, 4UL, MPFR_RNDN);
mpfr_log(tmps.mpfr, tmps.mpfr, MPFR_RNDN);
mpfr_add_ui(tmps.mpfr, tmps.mpfr, 1UL, MPFR_RNDN);
mpfr_atan(tmps.mpfr, tmps.mpfr, MPFR_RNDN);
mpfr_div(tmps.mpfr, tmps.mpfr, tmps.pi, MPFR_RNDN);
mpfr_add_d(tmps.mpfr, tmps.mpfr, 0.5, MPFR_RNDN);
// get brightness in steps of 1/256-ths
mpfr_mul_ui(tmps.mpfr, tmps.mpfr, 256, MPFR_RNDN);
uint8_t brightness = mpfr_get_ui(tmps.mpfr, MPFR_RNDZ);
*colour = (struct rgba_pix){ brightness, brightness, brightness, 0xff };
}
inline static struct rgba_pix interpolate(double t, struct rgba_pix a, struct rgba_pix b)
{
return (struct rgba_pix){
(double)a.r + t * ((double)b.r - (double)(b.r)),
(double)a.g + t * ((double)b.g - (double)(b.g)),
(double)a.b + t * ((double)b.b - (double)(b.b)),
0xff
};
}
void get_colour_cividis(struct rgba_pix *colour, const mpc_t val, struct temporaries tmps)
{
mpc_arg(tmps.mpfr, val, MPFR_RNDN);
mpfr_div(tmps.mpfr, tmps.mpfr, tmps.pi, MPFR_RNDN);
mpfr_mul_ui(tmps.mpfr, tmps.mpfr, 256, MPFR_RNDN);
uint8_t idx = mpfr_get_ui(tmps.mpfr, MPFR_RNDZ);
if (idx == 0xff) {
*colour = cividis_lut[idx];
return;
}
// interpolate linearly
mpfr_frac(tmps.mpfr, tmps.mpfr, MPFR_RNDN);
double t = mpfr_get_d(tmps.mpfr, MPFR_RNDN);
*colour = interpolate(t, cividis_lut[idx], cividis_lut[idx + 1]);
}
/* SPDX-License-Identifier: GPL-3.0-or-later */
#include "common.h"
struct temporaries init_temporaries(uint64_t prec)
{
struct temporaries tmps = { 0 };
mpfr_init2(tmps.pi, prec);
mpfr_const_pi(tmps.pi, MPFR_RNDN);
mpfr_init2(tmps.mpfr, prec);
mpc_init2(tmps.mpc, prec);
return tmps;
}
void destroy_temporaries(struct temporaries tmps)
{
mpfr_clear(tmps.pi);
mpfr_clear(tmps.mpfr);
mpc_clear(tmps.mpc);
}
/* SPDX-License-Identifier: GPL-3.0-or-later */
#include "eval.h"
void horner_eval(mpc_t res, const mpc_t z, const struct polynomial p)
{
mpc_set(res, p.coeffs[p.degree], MPFR_RNDN);
for (uint64_t i = p.degree; i > 0; i--) {
mpc_mul(res, res, z, MPFR_RNDN);
mpc_add(res, res, p.coeffs[i - 1], MPFR_RNDN);
}
}
/* SPDX-License-Identifier: GPL-3.0-or-later */
#include <stdbool.h>
#include <stdlib.h>
#include <ctype.h>
#include <errno.h>
#include <fcntl.h>
#include <unistd.h>
#include "common.h"
#include "colour.h"
#include "eval.h"
#include "png_simplified.h"
#include "config.h"
#define POLY_SIZE 20
struct options {
bool monochrome;
uint64_t prec;
uint64_t resolution;
mpc_t offset;
mpc_t x_step;
mpc_t y_step;
mpfr_t radius;
};
static inline void get_pos(const struct options opts, mpc_t z, uint64_t i, uint64_t j, struct temporaries tmps)
{
mpc_set(z, opts.offset, MPFR_RNDN);
mpc_mul_ui(tmps.mpc, opts.y_step, i, MPFR_RNDN);
mpc_add(z, z, tmps.mpc, MPFR_RNDN);
mpc_mul_ui(tmps.mpc, opts.x_step, j, MPFR_RNDN);
mpc_add(z, z, tmps.mpc, MPFR_RNDN);
}
struct options init_options(bool monochrome, uint64_t prec, uint64_t resolution)
{
struct options opts = { 0 };
opts.monochrome = monochrome;
opts.prec = prec;
opts.resolution = resolution;
mpc_init2(opts.offset, prec);
mpc_init2(opts.x_step, prec);
mpc_init2(opts.y_step, prec);
mpfr_init2(opts.radius, prec);
return opts;
}
void destroy_options(struct options opts)
{
mpc_clear(opts.offset);
mpc_clear(opts.x_step);
mpc_clear(opts.y_step);
mpfr_clear(opts.radius);
}
int get_uint(const char *s, uint64_t *n)
{
char *endptr;
int64_t n_signed;
n_signed = strtoll(s, &endptr, 10);
if (!*optarg || *endptr) {
return -1;
}
if ((n_signed == LLONG_MIN || n_signed == LLONG_MAX) && errno == ERANGE) {
return -1;
}
if (n_signed < 0) {
return -1;
}
*n = n_signed;
return 0;
}
int read_polynomial(struct options opts, const char *coefficient_file, struct polynomial *p)
{
int fd = open(coefficient_file, O_RDONLY);
if (fd < 0) {
perror("open");
return -1;
}
const uint64_t len = 1 << 12;
char *buf = calloc(len, sizeof(char));
if (!buf) {
perror("calloc");
return -1;
}
uint64_t degree = POLY_SIZE;
p->coeffs = calloc(degree, sizeof(mpc_t));
if (!p->coeffs) {
perror("calloc");
free(buf);
return -1;
}
off_t off = 0;
ssize_t n;
char *endptr;
uint64_t idx = 0;
while ((n = pread(fd, buf, len - 1, off)) > 0) {
if (n == 1 && isspace(buf[0])) {
break;
}
buf[n] = '\0';
mpc_init2(p->coeffs[idx], opts.prec);
if (mpc_strtoc(p->coeffs[idx], buf, &endptr, 10, MPFR_RNDN)) {
fprintf(stderr, "failed to read complex number from `%s'.\nMaybe the precission is too low.", buf);
for (uint64_t i = 0; i <= idx; i++) {
mpc_clear(p->coeffs[i]);
}
free(p->coeffs);
free(buf);
close(fd);
return -1;
}
off += endptr - buf;
idx++;
if (idx == degree) {
mpc_t *tmp = realloc(p->coeffs, (degree + POLY_SIZE) * sizeof(mpc_t));
if (!tmp) {
perror("realloc");
for (uint64_t i = 0; i < degree; i++) {
mpc_clear(p->coeffs[i]);
}
free(p->coeffs);
free(buf);
close(fd);
return -1;
}
p->coeffs = tmp;
degree += POLY_SIZE;
}
}
if (errno) {
perror("pread");
for (uint64_t i = 0; i < idx; i++) {
mpc_clear(p->coeffs[i]);
}
free(p->coeffs);
free(buf);
close(fd);
return -1;
}
if (idx == 0) {
fprintf(stderr, "coefficient file `%s' is empty\n", coefficient_file);
for (uint64_t i = 0; i < idx; i++) {
mpc_clear(p->coeffs[i]);
}
free(p->coeffs);
free(buf);
close(fd);
return -1;
}
p->degree = idx - 1;
free(buf);
close(fd);
return 0;
}
void destroy_polynomial(struct polynomial p)
{
for (uint64_t i = 0; i <= p.degree; i++) {
mpc_clear(p.coeffs[i]);
}
free(p.coeffs);
}
void fill(struct pic pic, struct options opts, const struct polynomial p)
{
void (*get_colour)(struct rgba_pix *, const mpc_t, struct temporaries) = opts.monochrome ? &get_colour_mono : &get_colour_cividis;
#pragma omp parallel
{
struct temporaries tmps = init_temporaries(opts.prec);
mpc_t z;
mpc_init2(z, opts.prec);
mpc_t val;
mpc_init2(val, opts.prec);
#pragma omp for
for (uint64_t i = 0; i < opts.resolution; i++) {
for (uint64_t j = 0; j < opts.resolution; j++) {
get_pos(opts, z, i, j, tmps);
mpc_abs(tmps.mpfr, z, MPFR_RNDN);
if (mpfr_greater_p(tmps.mpfr, opts.radius)) {
pic.pix[i * opts.resolution + j] = (struct rgba_pix){ 0, 0, 0, 0 };
continue;
}
horner_eval(val, z, p);
get_colour(pic.pix + i * opts.resolution + j, val, tmps);
}
}
destroy_temporaries(tmps);
mpc_clear(z);
mpc_clear(val);
mpfr_free_cache2(MPFR_FREE_LOCAL_CACHE);
}
}
void print_help(const char *exe)
{
printf(
"Usage:\n"
"\t%s [options] <coefficient_file>\n"
"\n"
"Options:\n"
" -h show this help message\n"
" -V show version information\n"
" -m use monochrome\n"
" -p PREC set precision\n"
" -c CENTER set center\n"
" -r RADIUS set radius\n"
" -R RESOLUTION set resolution\n"
" -o OUTFILE set output file\n",
exe);
}
void print_version(void)
{
printf(
"polyplot version %s\n"
"Copyright (C) 2022 Maik Herbers\n"
"\n"
"This program is free software: you can redistribute it and/or modify\n"
"it under the terms of the GNU General Public License as published by\n"
"the Free Software Foundation, either version 3 of the License, or\n"
"(at your option) any later version.\n"
"\n"
"This program is distributed in the hope that it will be useful,\n"
"but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
"MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
"GNU General Public License for more details.\n"
"\n"
"You should have received a copy of the GNU General Public License\n"
"along with this program. If not, see <https://www.gnu.org/licenses/>.\n",
POLYPLOT_VERSION);
}
int main(int argc, char **argv)
{
uint64_t prec = 53UL;
uint64_t resolution = 100UL;
const char *center_str = "(0 0)";
const char *radius_str = "1";
const char *out_file = "out.png";
const char *coefficient_file = "";
bool monochrome = false;
struct options opts;
while (true) {
switch (getopt(argc, argv, "+hVmp:c:r:R:o:")) {
case -1:
goto escape;
case 'h':
print_help(argv[0]);
return EXIT_SUCCESS;
case 'm':
monochrome = true;
break;
case 'p':
if (get_uint(optarg, &prec)) {
printf("invalid precision `%s'\n", optarg);
return EXIT_FAILURE;
}
if (prec < MPFR_PREC_MIN || prec > MPFR_PREC_MAX) {
printf("precision `%lu' outside of allowed range [%d, %ld]\n", prec, MPFR_PREC_MIN, MPFR_PREC_MAX);
return EXIT_FAILURE;
}
break;
case 'c':
center_str = optarg;
break;
case 'r':
radius_str = optarg;
break;
case 'R':
if (get_uint(optarg, &resolution)) {
printf("invalid resolution `%s'\n", optarg);
return EXIT_FAILURE;
}
break;
case 'o':
out_file = optarg;
break;
case 'V':
print_version();
return EXIT_SUCCESS;
case '?':
default:
print_help(argv[0]);
return EXIT_FAILURE;
}
}
escape:
if (argc - optind != 1) {
print_help(argv[0]);
return EXIT_FAILURE;
}
coefficient_file = argv[optind];
mpc_t tmp;
mpc_t center;
mpc_init2(tmp, prec);
mpc_init2(center, prec);
opts = init_options(monochrome, prec, resolution);
if (mpc_set_str(center, center_str, 10, MPFR_RNDN)) {
printf("invalid center argument `%s'\n", center_str);
return EXIT_FAILURE;
}
if (mpfr_set_str(opts.radius, radius_str, 10, MPFR_RNDN)) {
printf("invalid radius argument `%s'\n", radius_str);
return EXIT_FAILURE;
}
if (!mpfr_regular_p(opts.radius)) {
printf("invalid radius `%s'\n", center_str);
return EXIT_FAILURE;
}
mpc_set_fr_fr(tmp, opts.radius, opts.radius, MPFR_RNDN);
mpc_sub(opts.offset, center, tmp, MPFR_RNDN);
mpc_set_fr(opts.x_step, opts.radius, MPFR_RNDN);
mpc_div_ui(opts.x_step, opts.x_step, resolution, MPFR_RNDN);
mpc_mul_ui(opts.x_step, opts.x_step, 2, MPFR_RNDN);
mpc_mul_i(opts.y_step, opts.x_step, 1, MPFR_RNDN);
mpc_clear(tmp);
mpc_clear(center);
struct pic pic = init_pic(opts.resolution, opts.resolution);
if (!pic.pix) {
perror("calloc");
return EXIT_FAILURE;
}
struct polynomial p;
if (read_polynomial(opts, coefficient_file, &p)) {
return EXIT_FAILURE;
}
fill(pic, opts, p);
destroy_polynomial(p);
destroy_options(opts);
mpfr_free_cache();
save_to_png(pic, out_file);
free_pic(pic);
return EXIT_SUCCESS;
}
# SPDX-License-Identifier: GPL-3.0-or-later
src = [
'main.c',
'colour.c',
'common.c',
'png_simplified.c',
'eval.c',
]
executable(
meson.project_name(),
src,
include_directories : inc,
dependencies : deps,
install : true
)
/* SPDX-License-Identifier: GPL-3.0-or-later */
#include "png_simplified.h"
#include <stdlib.h>
struct pic init_pic(uint64_t height, uint64_t width)
{
png_image im = {
.version = PNG_IMAGE_VERSION,
.height = height,
.width = width,
.format = PNG_FORMAT_RGBA,
.opaque = NULL,
.flags = 0
};
struct pic pic = {
.im = im,
.pix = NULL
};
struct rgba_pix *pix = calloc(sizeof(struct rgba_pix), height * width);
if (!pix) {
return pic;
}
pic.pix = pix;
return pic;
}
void save_to_png(struct pic pic, const char *out)
{
png_image_write_to_file(&pic.im, out, 0, pic.pix, 0, NULL);
}
void free_pic(struct pic pic)
{
free(pic.pix);
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment