Miroirs/DB48X-on-DM42
1
0
Fork 0
mirror of https://github.com/c3d/DB48X-on-DM42.git synced 2024-09-28 03:20:53 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

polynomials: Add support for polynomials as a separate type

Browse source 
Enough algorithms need a separate notion of polynomials.
Devise a space-efficient representation for them.

Signed-off-by: Christophe de Dinechin <christophe@dinechin.org>
This commit is contained in:
Christophe de Dinechin 2024-04-25 22:15:55 +02:00
parent 86a2d5ae18
commit b05453ad25
6 changed files with 417 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
Makefile
View file

@ -260,6 +260,7 @@ CXX_SOURCES += \
src/menu.cc \
src/object.cc \
src/plot.cc \
src/polynomial.cc \
src/program.cc \
src/renderer.cc \
src/runtime.cc \

1
sim/db48x.pro
View file

@ -74,6 +74,7 @@ SOURCES += \
../src/menu.cc \
../src/object.cc \
../src/plot.cc \
../src/polynomial.cc \
../src/program.cc \
../src/renderer.cc \
../src/runtime.cc \

1
src/ids.tbl
View file

@ -1267,6 +1267,7 @@ ID(sparse_font)
ID(dmcp_font)
ID(tag)
ID(polynomial)
ID(conditional)
ID(while_conditional)

1
src/object.cc
View file

@ -60,6 +60,7 @@
#include "menu.h"
#include "parser.h"
#include "plot.h"
#include "polynomial.h"
#include "program.h"
#include "renderer.h"
#include "runtime.h"

241
src/polynomial.cc Normal file
View file

@ -0,0 +1,241 @@
// ****************************************************************************
// polynomial.c DB48X project
// ****************************************************************************
//
// File Description:
//
// Dense representation of multivariate polynomials
//
// Some operations on polynomials are much easier or faster if done
// with a numerical representation of the coefficients.
// We choose a dense representation here in line with the primary objective
// of DB48X to run on very memory-constrainted machines like the DM42
//
//
//
// ****************************************************************************
// (C) 2024 Christophe de Dinechin <christophe@dinechin.org>
// This software is licensed under the terms outlined in LICENSE.txt
// ****************************************************************************
// This file is part of DB48X.
//
// DB48X is free software: you can redistribute it and/or modify
// it under the terms outlined in the LICENSE.txt file
//
// DB48X is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// ****************************************************************************
#include "polynomial.h"
#include "arithmetic.h"
#include "expression.h"
#include "grob.h"
#include "integer.h"
size_t polynomial::variables() const
// ----------------------------------------------------------------------------
// Return the number of variables
// ----------------------------------------------------------------------------
{
byte_p p = byte_p(this);
byte_p first = p;
size_t length = leb128<size_t>(p);
size_t nvars = leb128<size_t>(p);
return (size_t(p - first) < length) ? nvars : 0;
}
symbol_g polynomial::variable(uint index) const
// ----------------------------------------------------------------------------
// Return the variable at the given index as a symbol
// ----------------------------------------------------------------------------
{
byte_p p = byte_p(this);
byte_p first = p;
size_t length = leb128<size_t>(p);
size_t nvars = leb128<size_t>(p);
if (index >= nvars)
return nullptr;
for (size_t v = 0; v < index; v++)
{
size_t vlen = leb128<size_t>(p);
p += vlen;
}
if (size_t(p - first) >= length)
return nullptr;
size_t vlen = leb128<size_t>(p);
return symbol::make(p, vlen);
}
PARSE_BODY(polynomial)
// ----------------------------------------------------------------------------
// No parsing for polynomials, they are only generated
// ----------------------------------------------------------------------------
{
return SKIP;
}
EVAL_BODY(polynomial)
// ----------------------------------------------------------------------------
// We can evaluate polynomials a bit faster than usual expressions
// ----------------------------------------------------------------------------
{
polynomial_g poly = o;
size_t nvars = poly->variables();
algebraic_g vars[nvars];
// Evaluate each of the variables exactly once (this is where we save time)
for (size_t v = 0; v < nvars; v++)
{
symbol_g var = poly->variable(v);
object_p evaluated = var->evaluate();
if (!evaluated)
return ERROR;
algebraic_g alg = evaluated->as_extended_algebraic();
if (!alg)
{
rt.type_error();
return ERROR;
}
vars[v] = alg;
}
// Loop over all factors
algebraic_g result = +integer::make(0);
for (auto term : *poly)
{
algebraic_g factor = term.factor();
for (size_t v = 0; v < nvars; v++)
{
ularge exponent = term.exponent();
algebraic_g value = +integer::make(exponent);
value = pow(vars[v], value);
factor = factor * value;
if (!factor)
return ERROR;
}
result = result + factor;
if (!result)
return ERROR;
}
// We are done, push the result
return rt.push(+result) ? OK : ERROR;
}
RENDER_BODY(polynomial)
// ----------------------------------------------------------------------------
// Render a polynomial as text
// ----------------------------------------------------------------------------
{
polynomial_g poly = o;
size_t nvars = poly->variables();
symbol_g vars[nvars];
// Get each of the variables
for (size_t v = 0; v < nvars; v++)
vars[v] = poly->variable(v);
// Loop over all factors
bool first = true;
unicode mul = Settings.UseDotForMultiplication() ? L'·' : L'×';
for (auto term : *poly)
{
// Separate terms with +
if (!first)
r.put('+');
first = false;
// Emit the factor
algebraic_g factor = term.factor();
factor->render(r);
for (size_t v = 0; v < nvars; v++)
{
ularge exponent = term.exponent();
if (exponent)
{
r.put(mul);
vars[v]->render(r);
if (exponent > 1)
{
r.put(unicode(L'↑'));
r.printf("%llu", exponent);
}
}
}
}
// We are done, push the result
return r.size();
}
GRAPH_BODY(polynomial)
// ----------------------------------------------------------------------------
// Render a polynomial as a graphic expression
// ----------------------------------------------------------------------------
{
polynomial_g poly = o;
size_t nvars = poly->variables();
grob_g vars[nvars];
// Get each of the variables and render it graphically
for (size_t v = 0; v < nvars; v++)
{
symbol_g sym = poly->variable(v);
grob_g var = sym->graph(g);
vars[v] = var;
}
// Loop over all factors
grob_g result = nullptr;
coord vr = 0;
cstring mul = Settings.UseDotForMultiplication() ? "·" : "×";
for (auto term : *poly)
{
// Render the factor
algebraic_g factor = term.factor();
grob_g factg = factor->graph(g);
coord vf = 0;
// Render the terms
for (size_t v = 0; v < nvars; v++)
{
ularge exponent = term.exponent();
if (exponent)
{
grob_g termg = vars[v];
coord vt = 0;
if (exponent > 1)
{
char exptxt[16];
snprintf(exptxt, sizeof(exptxt), "%llu", exponent);
termg = suscript(g, vt, termg, 0, exptxt);
vt = g.voffset;
}
factg = infix(g, vf, factg, 0, mul, vt, termg);
vf = g.voffset;
}
}
// Addition of terms
if (result)
result = infix(g, vr, result, 0, "+", vf, factg);
else
result = factg;
vr = g.voffset;
}
// We are done, push the result
return result;
}

172
src/polynomial.h Normal file
View file

@ -0,0 +1,172 @@
#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H
// ****************************************************************************
// polynomial.h DB48X project
// ****************************************************************************
//
// File Description:
//
// Dense representation of multivariate polynomials
//
// Some operations on polynomials are much easier or faster if done
// with a numerical representation of the coefficients.
// We choose a dense representation here in line with the primary objective
// of DB48X to run on very memory-constrainted machines like the DM42
//
//
//
// ****************************************************************************
// (C) 2024 Christophe de Dinechin <christophe@dinechin.org>
// This software is licensed under the terms outlined in LICENSE.txt
// ****************************************************************************
// This file is part of DB48X.
//
// DB48X is free software: you can redistribute it and/or modify
// it under the terms outlined in the LICENSE.txt file
//
// DB48X is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// ****************************************************************************
//
/* A polynomial is represented as a structure very similar to arrays,
except that the program block is replaced with an array block.
0. ID_polynomial
1. Total length for fast skipping
2. Number of variables
3. Sequence of variable names, each one being
3.1 Variable 1 name length
3.2 Variable 1 name
4. Length of array block
5. Sequence of array objects, each being in the form:
5.1. Factor value, a real or complex number
5.2. N variable exponents, one per variable
Exponents are sorted in decreasing lexicographic order.
For example 2/3 * X^37 * Y^42 + 1.25 * X^23 * Y^55 + (2+3i)*X - 1 is:
0. ID_polynomials
1. [Total length]
2. 2 (two variables, X and Y)
3. Two variables, Y and X (Y comes first because 55 > 37)
1 Y 1 X
4. [Array block length]
5. Decimal(1.25) 55 23
Fraction(2/3) 42 37
Complex(2+3i) 0 1
Neg_Integer(-1) 0 0
Polynomials are never parsed directly, but they can be built by symbolic
operations on expressions.
*/
#include "expression.h"
#include "symbol.h"
GCP(polynomial);
struct polynomial : expression
// ----------------------------------------------------------------------------
// Representation for polynomials
// ----------------------------------------------------------------------------
{
polynomial(id type, gcbytes bytes, size_t len)
: expression(type, bytes, len)
{ }
// Return total length of the polynomial in bytes
size_t length() const
{
return text::length();
}
// Access variables in the polynomial
size_t variables() const;
symbol_g variable(uint index) const;
// Iterating over factors and exponents
struct iterator
{
typedef iterator &value_type;
explicit iterator(polynomial_p poly, bool at_end = false)
: poly(poly), length(), variables(), offset()
{
byte_p p = byte_p(poly);
byte_p first = p;
length = leb128<size_t>(p);
variables = leb128<size_t>(p);
for (size_t v = 0; v < variables; v++)
{
// Skip each name
size_t vlen = leb128<size_t>(p);
p += vlen;
}
offset = at_end ? length : p - first;
}
algebraic_p factor()
{
algebraic_p scalar = algebraic_p(poly) + offset;
object_p exponents = scalar->skip();
offset = exponents - object_p(poly);
return scalar;
}
ularge exponent()
{
byte_p p = byte_p(poly) + offset;
uint exp = leb128<ularge>(p);
offset = p - byte_p(poly);
return exp;
}
bool operator==(const iterator &o) const
{
return +o.poly == +poly
&& o.offset == offset
&& o.length == length
&& o.variables == variables;
}
bool operator!=(const iterator &o) const
{
return !(o==*this);
}
// We don't increment, because it's really factor() and exponent()
// that increment the offset.
iterator& operator++()
{
return *this;
}
iterator operator++(int)
{
return *this;
}
value_type operator*()
{
return *this;
}
polynomial_g poly;
size_t length;
size_t variables;
size_t offset;
};
iterator begin() const { return iterator(this); }
iterator end() const { return iterator(this, true); }
public:
OBJECT_DECL(polynomial);
PARSE_DECL(polynomial);
EVAL_DECL(polynomial);
RENDER_DECL(polynomial);
GRAPH_DECL(polynomial);
};
#endif // POLYNOMIAL_H