Skip to content

Instantly share code, notes, and snippets.

@ochafik
Last active February 23, 2022 01:03
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ochafik/8801e5806a362ab3037ff423289f5348 to your computer and use it in GitHub Desktop.
Save ochafik/8801e5806a362ab3037ff423289f5348 to your computer and use it in GitHub Desktop.
WIP CGAL Cached Number for Kernels

NOT READY TO POST YET This is still speculative / not backed by concrete results yet.

Trading space for speed in exact computations

Complex models can take time to render in OpenSCAD, even with the 100x speedups achievable with the experimental fast-csg option.

Lots of that rendering time is spent doing expensive computations on exact numbers. Within the same render, common subparts are cached, but between runs, everything starts from scratch.

There are two possible leads to cache computations:

  1. Persist the rendering cache across between runs
  2. Cache the results at the number level

While 1. is relatively straightforward (e.g. should be able to dump geometry in a binary format to disk), 2. is still highly speculative and the subject of this post.

Singleton numbers?

The basic idea is to only ever allocate a single instance of the "same" number and cache the results of all its operations with all other numbers (and since we use exact numeric types in OpenSCAD, same means really really exactly the same). Said results being themselves pointers to singleton numbers (or booleans).

This isn't really an option in general numerical processing, as there are way too many double values to cache, but in the case of OpenSCAD there's only so many numbers floating around, if I may say (see next section).

We currently use two different kernels

  • Cartesian<Gmpq>: our legacy basic no-thrills Gmpq kernel. Does every operation it's asked to do, with exact accuracy.
  • Epeck (used by fast-csg along with fast corefinement operations): a lazy filtering kernel, keeps a tree of operations that does a lot of voodoo, avoiding actual operations in many cases (using "filtering" to do double operations on bounds and defer actual operations to later - sometimes never - hence the name). We actually have to undo some of that voodoo in the fast-csg code to force numbers to not be too lazy, otherwise performance can degrade dramatically.

While Epeck limits the number of operations needed and fast-tracks some, I propose a different (potentially complementary) approach: do all the operations needed, but only once, ever. And canonicalize numbers and all the potential results of all of their operations with all potential other numbers.

(we can combine those canonicalized numbers with Epeck-style filtering, of course, to do less actual operations on these singletons)

A matter of scale: how many numbers are there?

Since OpenSCAD isn't running in an infinite platonic space of abstract numbers and has to deal with the vicissitudes of finite memory, let's put a price tag on this approach.

Conceptually, you might think the number of unique numbers floating around is O(V), where V is the number of vertices that were ever created during the rendering. This is true with original numbers (say, the ones in the STL meshes imported), and with numbers created during operations like transforms or extrusions, etc., but then CSG operations can create non linear amounts of by-products.

But the good news is: the storage complexity of cached singleton operations with each other should be the same as the current execution complexity (if your model renders in quadratic time, you'll need quadratic growth of cache / memory storage), and it will reduces execution time by a massive constant (leaving its complexity unchanged), swapping expensive exact number computations with integer arithmetics (hashtable lookups).

With examples/Basics/CSG.html here are the stats (nef mode, no fast-csg, typedef CGAL::Cartesian<SingletonNumber<CGAL::Gmpq>> CGAL_Kernel3: (note: 0 and 1 are fast-tracked but would otherwise be by far the most resolved values)

values_: 63731
doubleValues_: 63731
idsByDoubleValue_: 63731
makeFractionOpsCache: 0
unaryMinusOpsCache: 3362
sqrtOpsCache: 0
plusOpsCache: 8383
minusOpsCache: 27028
timesOpsCache: 77001
divideOpsCache: 497
minOpsCache: 0
maxOpsCache: 0
equalOpsCache: 1027
lessOpsCache: 19970
greaterOpsCache: 3277
Common resolved doubles (60211):
  0: 	resolved 2051 times.
  2: 	resolved 204 times.
  -2.07912: 	resolved 105 times.
  2.07912: 	resolved 105 times.
  -4.06737: 	resolved 104 times.
  4.06737: 	resolved 104 times.
  -5.87785: 	resolved 100 times.
  5.87785: 	resolved 100 times.
  -9.94522: 	resolved 97 times.
  9.94522: 	resolved 97 times.
  -9.51057: 	resolved 96 times.
  -8.66025: 	resolved 96 times.
  -7.43145: 	resolved 96 times.
  7.43145: 	resolved 96 times.
  8.66025: 	resolved 96 times.
  9.51057: 	resolved 96 times.
  -7.5: 	resolved 91 times.
  7.5: 	resolved 90 times.
  2.03997: 	resolved 57 times.
  -0.131264: 	resolved 54 times.

With the example corefinement snippet below, getting the following output (note: 0 is fast-tracked but would otherwise be by far the most resolved value)

corefinement: true
values_: 2463
doubleValues_: 2463
idsByDoubleValue_: 2463
makeFractionOpsCache: 0
unaryMinusOpsCache: 351
sqrtOpsCache: 0
plusOpsCache: 641
minusOpsCache: 714
timesOpsCache: 2464
divideOpsCache: 29
minOpsCache: 0
maxOpsCache: 0
equalOpsCache: 82
lessOpsCache: 853
greaterOpsCache: 405
Common resolved doubles (2463):
  1: 	resolved 488 times.
  0: 	resolved 102 times.
  -0.382683: 	resolved 22 times.
  0.382683: 	resolved 18 times.
  -0.92388: 	resolved 16 times.
  0.92388: 	resolved 16 times.
  -1: 	resolved 14 times.
  -0.353554: 	resolved 14 times.
  -0.146447: 	resolved 14 times.
  0.270598: 	resolved 14 times.
  -0.270598: 	resolved 13 times.
  0.146447: 	resolved 13 times.
  0.353554: 	resolved 13 times.
  -0.653281: 	resolved 12 times.
  0.0606601: 	resolved 11 times.
  0.146446: 	resolved 11 times.
  0.146447: 	resolved 11 times.
  0.653281: 	resolved 11 times.
  504: 	resolved 11 times.
  -0.146447: 	resolved 10 times.

A concrete example (TODO)

A very complex rendering of a heavy model in OpenSCAD takes X minutes and eats up a peak Y MB of RAM. It creates Z instances of Gmpq numbers (in fast-csg mode, using the filtered Epeck kernel), which can be canonicalized to W unique instances. It does M operations

With SingletonNumber<Gmpq>, we only ever create W instances and cache N operation results. This takes Y' MB of RAM. And we only do M' operations, which results in a run time of T' minutes (?% improvement).

Bring OpenSCAD fully to exact world

Currently we use doubles in many places (look for usages of Epick and PolySet in the codebase), either for speed or backwards compatibility:

  • Transforms
  • Hull
  • Minkowski
  • (maybe more 2D operations?)

We'd need to use exact numbers through out to make the best use of the singleton numbers. Instead of numerical operations, the state of the program can be seen as transitioning on a large-scale numeric transition graph, with starting and ending points being the double <-> exact conversions.

Saving electricity?

With a hot persisted cache, even a large rendering wouldn't need a single Gmpq (or double) operation!!

For each heavy model, one could precompute and share such a cache fiel on/with their rendering farm to cut computation costs, save electricity and time. No need to compute 1 + 1 with exact numbers anymore, ever.

This only holds true if looking up the results of 1 + 1 is cheaper than the exact computation itself, of course.

How complex are Gmpq operations

Operation on exact rationals (see internals) involve separate operations on numerator and denominator (themselves variable length big integers) which are moderately slow, followed by GCD operations (likely the slowest bit) to ensure the final numerator and denominator have no common divider. Creating results for each big integer involves allocating memory (for a variable amount of limbs in GMP parlance).

In contrast, operations on "hot" singletons are just fixed-time hash operations and memory access (hopefully cache-aided).

MRU Map and cache locality

If we're gonna keep numbers around forever, we might run out of space and want to ditch the unused numbers at some point.

At the beginning of each render, we increment an runId number. Each time we cache-hit a number, we 'touch' it by updating its lastRunIdRead. At the end of the render, we can ditch the older entries based on that number.

For bonus points we can compact the array of values to put freshly used values (their double conversion mostly) close to each other to improve cache locality of the final exports.

We can serialize the numbers cache to disk, maybe even use some memory mapping to have the hashtable file live in memory with no extra cost (a file for exact Gmpq values, one for the double values, one for the hashtable).

Implementation notes

  • Singleton property not absolutely guaranteed unless two equal numbers are guaranteed to convert to the same double value with CGAL::to_double. If they don't, we'll get some cache misses. Moderate amounts of such occurrences aren't bothering, but need to check if this can happen at all / how often.
  • NumberId: is a single non-zero size_t that indexes both exact values and double conversions
/*
g++ -std=c++14 \
-DDEBUG=1 -O0 -g \
-DCGAL_USE_GMPXX=1 -lgmp -lmpfr \
-o sing src/cgal-singleton-number.cc
*/
#include <CGAL/Gmpq.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Nef_polyhedron_3.h>
#include <CGAL/Cartesian.h>
#include "cgal-singleton-number.h"
#include <fstream>
extern std::string example1; // see bottom of file
extern std::string example2;
// template <class T>
// class Foo {
// public:
// static int bar;
// int getBar() {
// return bar;
// }
// };
// template<>
// int Foo<double>::bar = 10;
// template <class UnderlyingType>
// class SingletonCache2 {
// public:
// static Foo<UnderlyingType> foo;
// // static SingletonCache2<UnderlyingType> instance;
// int getBar() {
// return Foo<UnderlyingType>::bar;
// // return foo.getBar();
// }
// };
// template <>
// Foo<double> SingletonCache2<double>::foo;
// template<typename FT>
// SingletonCache<FT> SingletonNumber<FT>::cache;
// template<>
// SingletonCache<CGAL::Gmpq> SingletonNumber<CGAL::Gmpq>::cache;
// template <>
// SingletonCache<CGAL::Gmpq> SingletonNumber<CGAL::Gmpq>::cache;
// template <typename UnderlyingType>
// SingletonCache<UnderlyingType>& getCache() {
// static SingletonCache<UnderlyingType> cache;
// return cache;
// }
// template <>
// SingletonCache<CGAL::Gmpq>& getCache<CGAL::Gmpq>();
// static SingletonCache<UnderlyingType> cache;
// template<typename T>
// class A
// {
// public:
// static T val;
// static void init()
// {
// val=0;
// }
// };
// template<typename T> T A<T>::val; // only change here
// int main()
// {
// // A::init();
// A<double>::init();
// return 0;
// }
typedef SingletonNumber<CGAL::Gmpq> FT;
typedef CGAL::Cartesian<FT> K;
// typedef CGAL::Cartesian<CGAL::Gmpq> K;
// typedef CGAL::Simple_cartesian<double> K;
// typedef CGAL::Epeck K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef CGAL::Polyhedron_3<K> Poly;
typedef CGAL::Nef_polyhedron_3<K> Nef;
// int main0()
// {
// typedef Mesh::Vertex_index vertex_descriptor;
// typedef Mesh::Face_index face_descriptor;
// Mesh m;
// // Add the points as vertices
// vertex_descriptor u = m.add_vertex(K::Point_3(0,1,0));
// vertex_descriptor v = m.add_vertex(K::Point_3(0,0,0));
// vertex_descriptor w = m.add_vertex(K::Point_3(1,1,0));
// vertex_descriptor x = m.add_vertex(K::Point_3(1,0,0));
// m.add_face(u,v,w);
// face_descriptor f = m.add_face(u,v,x);
// if(f == Mesh::null_face())
// {
// std::cerr<<"The face could not be added because of an orientation error."<<std::endl;
// f = m.add_face(u,x,v);
// assert(f != Mesh::null_face());
// }
// SingletonNumber<CGAL::Gmpq>::cache.debug();
// return 0;
// }
int main()
{
try {
#ifdef TEST
bool b;
FT x(1), y(2), z;
// z = x + y;
// z = x - y;
// z = x * y;
// z = x / y;
// z += x;
// z -= x;
// z *= x;
// z /= x;
// z = -x;
// b = x < y;
// b = x <= y;
// b = x > y;
// b = x >= y;
// z = 2 + x;
// z = 2 - x;
// z = 2 * x;
// z = 2 / x;
// z = x + 2;
// z = x - 2;
// z = x * 2;
// z = x / 2;
auto expect = [&](auto actual, auto expected) {
if ((actual != expected) || !(actual == expected)) {
std::cout << "actual == " << CGAL::to_double(actual) << " (!= expected " << expected << ")\n";
assert(z == expected);
}
};
auto expectBool = [&](bool actual, bool expected) {
if (actual != expected) {
std::cout << "actual == " << (actual ? "true" : "false") << " (!= expected " << (expected ? "true" : "false") << ")\n";
assert(z == expected);
}
};
//
// Plus
//
expect(FT(0) + 1, 1);
expect(2 + FT(0), 2);
expect(FT(10) + 0, 10);
expect(0 + FT(11), 11);
expect(FT(2) + 3, 5);
expect(2 + FT(4), 6);
expect(FT(2) + FT(5), 7);
//
// Minus
//
expect(FT(0) - 1, -1);
expect(2 - FT(0), 2);
expect(FT(20) - 0, 20);
expect(0 - FT(20), -20);
expect(FT(2) + 3, 5);
expect(2 + FT(4), 6);
expect(FT(2) + FT(5), 7);
// Zero equality
expect(FT(0), 0);
expect(FT(-0), 0);
expect(FT(0), FT(-0));
expect(FT(1) - FT(1), 0);
expect(FT(0) - FT(1), FT(-1));
//
// Multiply
//
expect(FT(0) * 2, 0);
expect(2 * FT(0), 0);
expect(0 * FT(2), 0);
expect(FT(2) * 0, 0);
expect(FT(2) * 3, 6);
expect(2 * FT(4), 8);
expect(FT(2) * FT(5), 10);
//
// Divide
//
expect(FT(0) / 2, 0);
expect(0 / FT(2), 0);
expect(0 / FT(0), 0);
expect(4 / FT(2), 2);
expect(FT(9) / 3, 3);
expect(FT(12) / FT(3), 4);
// CGAL::Gmpq(1, 0) / CGAL::Gmpq(0, 1);
// FT(1) / FT(0);
//
// Less
//
expectBool(FT(0) < 2, true);
expectBool(FT(0) < -1, false);
expectBool(FT(0) < 0, false);
expectBool(FT(2) < 3, true);
expectBool(FT(2) < 0, false);
expectBool(FT(-1) < 0, true);
expectBool(FT(3) < 3, false);
expectBool(FT(10) < FT(11), true);
expectBool(FT(11) < FT(10), false);
z = -y;
expect(z, -2);
z = 3.0 + y;
expect(z, 5);
z = 3.0 - y;
expect(z, 1);
// z = 3.0 * y;
// expect(z, 6);
// z = 3.0 / y;
// expect(z, CGAL::Gmpq(3, 2));
z = x + 3.0;
expect(z, 4);
z = x - 3.0;
expect(z, -2);
// z = x * 3.0;
// expectZ(z, 3);
// z = x / 3.0;
// expectZ(z, CGAL::Gmpq(1, 3));
// std::cout << SingletonNumber<CGAL::Gmpq>::cache.values_.size() << " values\n";
// for (auto &value : SingletonNumber<CGAL::Gmpq>::cache.values_) {
// std::cout << value << "\n";
// }
#else
// Nef nef1, nef2;
Poly poly1, poly2;
Mesh mesh1, mesh2;
{
std::stringstream ss;
ss << example1;
ss >> poly1;
ss >> mesh1;
assert(CGAL::is_closed(mesh1)); // OK
assert(CGAL::is_valid_halfedge_graph(mesh1, /* verb */ false)); // OK
}
{
std::stringstream ss;
ss << example2;
ss >> poly2;
ss >> mesh2;
assert(CGAL::is_closed(mesh2)); // OK
assert(CGAL::is_valid_halfedge_graph(mesh2, /* verb */ false)); // OK
}
CGAL::Polygon_mesh_processing::triangulate_faces(mesh1);
CGAL::Polygon_mesh_processing::triangulate_faces(mesh2);
std::cout << "corefinement: " << (CGAL::Polygon_mesh_processing::corefine_and_compute_intersection(mesh1, mesh2, mesh1) ? "true" : "false") << "\n";
Nef nef1(poly1), nef2(poly2);
nef1 += nef2;
#endif
} catch (CGAL::Failure_exception& e) {
std::cerr << "CGAL error: " << e.what() << "\n";
}
SingletonNumber<CGAL::Gmpq>::cache.debug();
return 0;
}
// `sphere(1, $fn=8);` in OpenSCAD
std::string example1(
"OFF\n\
32 60 0\n\
0.270598 0.270598 0.92388\n\
0 0.382683 0.92388\n\
-0.270598 0.270598 0.92388\n\
-0.382683 0 0.92388\n\
-0.270598 -0.270598 0.92388\n\
0 -0.382683 0.92388\n\
0.270598 -0.270598 0.92388\n\
0.382683 0 0.92388\n\
0.653281 0.653281 0.382683\n\
0 0.92388 0.382683\n\
-0.653281 0.653281 0.382683\n\
-0.92388 0 0.382683\n\
-0.653281 -0.653281 0.382683\n\
0 -0.92388 0.382683\n\
0.653281 -0.653281 0.382683\n\
0.92388 0 0.382683\n\
0.653281 0.653281 -0.382683\n\
0 0.92388 -0.382683\n\
-0.653281 0.653281 -0.382683\n\
-0.92388 0 -0.382683\n\
-0.653281 -0.653281 -0.382683\n\
0 -0.92388 -0.382683\n\
0.653281 -0.653281 -0.382683\n\
0.382683 0 -0.92388\n\
0.92388 0 -0.382683\n\
0.270598 0.270598 -0.92388\n\
0 0.382683 -0.92388\n\
-0.270598 0.270598 -0.92388\n\
-0.382683 0 -0.92388\n\
-0.270598 -0.270598 -0.92388\n\
0 -0.382683 -0.92388\n\
0.270598 -0.270598 -0.92388\n\
3 7 15 8\n\
3 7 8 0\n\
3 0 8 9\n\
3 0 9 1\n\
3 1 9 10\n\
3 1 10 2\n\
3 2 10 11\n\
3 2 11 3\n\
3 3 11 12\n\
3 3 12 4\n\
3 4 12 13\n\
3 4 13 5\n\
3 5 13 14\n\
3 5 14 6\n\
3 6 14 15\n\
3 6 15 7\n\
3 15 24 16\n\
3 15 16 8\n\
3 8 16 17\n\
3 8 17 9\n\
3 9 17 18\n\
3 9 18 10\n\
3 10 18 19\n\
3 10 19 11\n\
3 11 19 20\n\
3 11 20 12\n\
3 12 20 21\n\
3 12 21 13\n\
3 13 21 22\n\
3 13 22 14\n\
3 14 22 24\n\
3 14 24 15\n\
3 24 23 25\n\
3 24 25 16\n\
3 16 25 26\n\
3 16 26 17\n\
3 17 26 27\n\
3 17 27 18\n\
3 18 27 28\n\
3 18 28 19\n\
3 19 28 29\n\
3 19 29 20\n\
3 20 29 30\n\
3 20 30 21\n\
3 21 30 31\n\
3 21 31 22\n\
3 22 31 23\n\
3 22 23 24\n\
3 0 1 7\n\
3 5 1 3\n\
3 3 1 2\n\
3 5 3 4\n\
3 7 1 5\n\
3 7 5 6\n\
3 30 29 28\n\
3 26 30 28\n\
3 27 26 28\n\
3 25 23 26\n\
3 23 30 26\n\
3 31 30 23\n\
");
// Derived from `sphere(1, $fn=3);` in OpenSCAD (edited to have integer coordinates, but same issue as the original + same general shape)
std::string example2(
"OFF\n\
6 8 0\n\
-3.0 6.0 7.0\n\
-3.0 -6.0 7.0\n\
7.0 0 -7.0\n\
7.0 0 7.0\n\
-3.0 6.0 -7.0\n\
-3.0 -6.0 -7.0\n\
3 1 3 0\n\
3 3 2 4\n\
3 3 4 0\n\
3 0 4 5\n\
3 0 5 1\n\
3 1 5 2\n\
3 1 2 3\n\
3 2 5 4\n\
");
// Portions of this file are Copyright 2021 Google LLC, and licensed under GPL2+. See COPYING.
//
// Based on my related "fuzzy number" attempt posted here: https://github.com/CGAL/cgal/issues/5490
//
#pragma once
#include <cmath>
#include <csignal>
#include <iostream>
#include <type_traits>
#include <unordered_map>
template <class FT>
class SingletonNumber;
template <class FT>
class SingletonCache {
typedef size_t NumberId;
typedef std::pair<NumberId, NumberId> NumberPair;
struct NumberPairHash
{
std::size_t operator() (const NumberPair &pair) const {
return std::hash<NumberId>()(pair.first) ^ std::hash<NumberId>()(pair.second);
}
};
typedef std::unordered_map<NumberId, NumberId> UnaryNumericOpCache;
typedef std::unordered_map<NumberPair, NumberId, NumberPairHash> BinaryNumericOpCache;
typedef std::unordered_map<NumberPair, bool, NumberPairHash> BinaryBooleanOpCache;
friend class SingletonNumber<FT>;
std::vector<FT> values_;
std::vector<double> doubleValues_;
std::unordered_multimap<double, NumberId> idsByDoubleValue_;
#ifdef DEBUG
std::unordered_map<double, size_t> debugDoubleResolutions_;
#endif
public:
SingletonCache() : values_(2), doubleValues_(2) {
// id=0 has value zero built by FT()
idsByDoubleValue_.insert(std::make_pair(0.0, 0));
idsByDoubleValue_.insert(std::make_pair(0.0, 1));
assert(values_[0] == FT(0.0));
assert(doubleValues_[0] == 0.0);
values_[1] = FT(1.0);
doubleValues_[1] = 1.0;
}
// SingletonCache<Gmpq>::makeFractionOpsCache takes 2 ids of Gmpz as key and gives a Gmpq value valid in SingletonCache<Gmpq>
BinaryNumericOpCache makeFractionOpsCache;
UnaryNumericOpCache unaryMinusOpsCache;
UnaryNumericOpCache sqrtOpsCache;
BinaryNumericOpCache plusOpsCache;
BinaryNumericOpCache minusOpsCache;
BinaryNumericOpCache timesOpsCache;
BinaryNumericOpCache divideOpsCache;
BinaryNumericOpCache minOpsCache;
BinaryNumericOpCache maxOpsCache;
BinaryBooleanOpCache equalOpsCache;
BinaryBooleanOpCache lessOpsCache;
BinaryBooleanOpCache greaterOpsCache;
void printStats() {
std::cout << "values_: " << values_.size() << "\n";
std::cout << "doubleValues_: " << doubleValues_.size() << "\n";
std::cout << "idsByDoubleValue_: " << idsByDoubleValue_.size() << "\n";
std::cout << "makeFractionOpsCache: " << makeFractionOpsCache.size() << "\n";
std::cout << "unaryMinusOpsCache: " << unaryMinusOpsCache.size() << "\n";
std::cout << "sqrtOpsCache: " << sqrtOpsCache.size() << "\n";
std::cout << "plusOpsCache: " << plusOpsCache.size() << "\n";
std::cout << "minusOpsCache: " << minusOpsCache.size() << "\n";
std::cout << "timesOpsCache: " << timesOpsCache.size() << "\n";
std::cout << "divideOpsCache: " << divideOpsCache.size() << "\n";
std::cout << "minOpsCache: " << minOpsCache.size() << "\n";
std::cout << "maxOpsCache: " << maxOpsCache.size() << "\n";
std::cout << "equalOpsCache: " << equalOpsCache.size() << "\n";
std::cout << "lessOpsCache: " << lessOpsCache.size() << "\n";
std::cout << "greaterOpsCache: " << greaterOpsCache.size() << "\n";
#ifdef DEBUG
std::vector<std::pair<double, size_t>> pairs(debugDoubleResolutions_.begin(), debugDoubleResolutions_.end());
struct {
bool operator()(const std::pair<double, size_t>& a, const std::pair<double, size_t>& b) const {
if (a.second == b.second) {
return a.first < b.first;
}
return a.second > b.second;
}
} order;
std::sort(pairs.begin(), pairs.end(), order);
std::cerr << "Common resolved doubles (" << pairs.size() << "):\n";
for (auto i = 0; i < pairs.size() && i < 20; i++) {
std::cerr << " " << pairs[i].first << ": \tresolved " << pairs[i].second << " times.\n";
}
#endif
}
NumberId getId(const FT& value) {
auto doubleValue = CGAL::to_double(value);
#ifdef DEBUG
debugDoubleResolutions_[doubleValue]++;
#endif
if (doubleValue == 0 && value == 0) {
// Fast track, no hashtable lookups.
return 0;
}
#if 1
auto it = idsByDoubleValue_.find(doubleValue);
while (it != idsByDoubleValue_.end()) {
if (it->first != doubleValue) {
break;
// throw 0;
}
// Multiple numbers can have the same double approximation but be actually different in exact precision.
// Exact (potentially very expensive) equality here.
auto existingId = it->second;
auto &existingValue = getValue(existingId);
if (existingValue == value) {
#ifdef DEBUG
if (CGAL::to_double(existingValue) != doubleValue) {
throw 0;
}
#endif
return existingId;
}
++it;
}
#endif
auto newId = values_.size();
idsByDoubleValue_.insert(it, std::make_pair(doubleValue, newId));
values_.push_back(value);
doubleValues_.push_back(doubleValue);
return newId;
}
const FT& getValue(NumberId id) {
return values_[id];
}
double getDoubleValue(NumberId id) {
return doubleValues_[id];
}
};
// template <typename FT>
// SingletonCache<FT>& getCache();
template <class FT>
class SingletonNumber {
typedef size_t NumberId;
typedef SingletonNumber<FT> Type;
// static SingletonCache<FT> cache;
NumberId id_;
// THe bool is just here to avoid any implicit usage & conflict with other ctors.
SingletonNumber(NumberId id, bool) : id_(id) {}
public:
static SingletonCache<FT> cache;
template <typename T, std::enable_if_t<std::is_arithmetic<T>::value, bool> = true>
SingletonNumber(const T& value) : id_(value == 0 ? 0 : value == 1 ? 1 : cache.getId(FT(value))) {} // static_cast<double>(value)?
SingletonNumber(double value) : id_(value == 0 ? 0 : value == 1 ? 1 : cache.getId(FT(value))) {}
SingletonNumber(const FT& value) : id_(cache.getId(value)) {}
SingletonNumber(const Type& other) : id_(other.id_) {}
SingletonNumber() : id_(0) {} //cache.getId(FT(0.0))) {}
const FT& underlyingValue() const {
return cache.getValue(id_);
}
explicit operator double() const {
return cache.getDoubleValue(id_);
}
bool isZero() const {
return id_ == 0;
}
bool isOne() const {
return id_ == 1;
}
Type &operator=(const Type& other) {
id_ = other.id_;
return *this;
}
#define SINGLETON_NUMBER_OP_TEMPLATE(T) \
template <typename T, std::enable_if_t<(std::is_arithmetic<T>::value || std::is_assignable<FT, T>::value), bool> = true>
template <class Operation, class CacheType>
static auto getCachedUnaryOp(NumberId operandId, CacheType &cache, Operation op) {
auto it = cache.find(operandId);
if (it != cache.end()) {
return it->second;
}
auto result = op(operandId);
cache.insert(it, std::make_pair(operandId, result));
return result;
}
template <class Operation, class CacheType>
static auto getCachedBinaryOp(NumberId lhsId, NumberId rhsId, CacheType &cache, Operation op) {
auto key = std::make_pair(lhsId, rhsId);
auto it = cache.find(key);
if (it != cache.end()) {
return it->second;
}
auto result = op(lhsId, rhsId);
cache.insert(it, std::make_pair(key, result));
return result;
}
template <class RawOperation>
struct UnaryNumericOp {
NumberId operator()(NumberId operand) {
RawOperation op;
return cache.getId(op(cache.getValue(operand)));
}
};
template <class RawOperation>
struct BinaryNumericOp {
NumberId operator()(NumberId lhs, NumberId rhs) {
RawOperation op;
return cache.getId(op(cache.getValue(lhs), cache.getValue(rhs)));
}
};
template <class RawOperation>
struct BinaryBooleanOp {
bool operator()(NumberId lhs, NumberId rhs) {
RawOperation op;
return op(cache.getValue(lhs), cache.getValue(rhs));
}
};
Type operator-() const {
struct Op {
FT operator()(const FT &operand) {
return -operand;
}
};
return Type(getCachedUnaryOp(id_, cache.unaryMinusOpsCache, UnaryNumericOp<Op>()), false);
}
Type sqrt() const {
struct Op {
FT operator()(const FT &operand) {
return std::sqrt(operand);
}
};
return Type(getCachedUnaryOp(id_, cache.sqrtOpsCache, UnaryNumericOp<Op>()), false);
}
bool operator==(const Type& other) const {
if (id_ == other.id_) {
return true;
}
if (isZero() != other.isZero()) {
return false;
}
if (isOne() != other.isOne()) {
return false;
}
struct Op {
bool operator()(const FT &lhs, const FT& rhs) {
return lhs == rhs;
}
};
return getCachedBinaryOp(id_, other.id_, cache.equalOpsCache, BinaryBooleanOp<Op>());
}
SINGLETON_NUMBER_OP_TEMPLATE(T) bool operator==(const T& x) const {
// Fast track to avoid resolving x's id.
if (x == 0) {
return isZero();
} else if (isZero()) {
return false;
}
if (x == 1) {
return isOne();
} else if (isOne()) {
return false;
}
return *this == Type(x);
}
bool operator<(const Type& other) const {
if (id_ == other.id_) {
return false;
}
struct Op {
bool operator()(const FT &lhs, const FT& rhs) {
return lhs < rhs;
}
};
return getCachedBinaryOp(id_, other.id_, cache.lessOpsCache, BinaryBooleanOp<Op>());
}
SINGLETON_NUMBER_OP_TEMPLATE(T) bool operator<(const T& x) const {
// Fast track to avoid resolving x's id.
if (x == 0 && isZero()) {
return false;
}
// if (x == 1 && isOne()) {
// return false;
// }
return *this < Type(x);
}
// a >= b if !(a < b)
bool operator>=(const Type& other) const {
return !(*this < other);
}
SINGLETON_NUMBER_OP_TEMPLATE(T) bool operator>=(const T& x) const {
return !(*this < x);
}
bool operator>(const Type& other) const {
if (id_ == other.id_) {
return false;
}
struct Op {
bool operator()(const FT &lhs, const FT& rhs) {
return lhs > rhs;
}
};
return getCachedBinaryOp(id_, other.id_, cache.greaterOpsCache, BinaryBooleanOp<Op>());
}
SINGLETON_NUMBER_OP_TEMPLATE(T) bool operator>(const T& x) const {
// Fast track to avoid resolving x's id.
if (x == 0 && isZero()) {
return false;
}
// if (x == 1 && isOne()) {
// return false;
// }
return *this > Type(x);
}
// a <= b if !(a > b)
bool operator<=(const Type& other) const {
return !(*this > other);
}
SINGLETON_NUMBER_OP_TEMPLATE(T) bool operator<=(const T& x) const {
return !(*this > x);
}
bool operator!=(const Type& other) const {
return !(*this == other);
}
SINGLETON_NUMBER_OP_TEMPLATE(T) bool operator!=(const T& x) const {
// Fast track to avoid resolving x's id.
if (x == 0) {
return !isZero();
} else if (isZero()) {
return false;
}
if (x == 1) {
return !isOne();
} else if (isOne()) {
return false;
}
return !(*this == Type(x));
}
Type min(const Type& other) const {
struct Op {
FT operator()(const FT &lhs, const FT& rhs) {
return std::min(lhs, rhs);
}
};
return Type(getCachedBinaryOp(id_, other.id_, cache.minOpsCache, BinaryNumericOp<Op>()), false);
}
Type max(const Type& other) const {
struct Op {
FT operator()(const FT &lhs, const FT& rhs) {
return std::max(lhs, rhs);
}
};
return Type(getCachedBinaryOp(id_, other.id_, cache.maxOpsCache, BinaryNumericOp<Op>()), false);
}
Type operator+(const Type& other) const {
if (isZero()) {
return other;
}
if (other.isZero()) {
return *this;
}
struct Op {
FT operator()(const FT &lhs, const FT& rhs) {
return lhs + rhs;
}
};
return Type(getCachedBinaryOp(id_, other.id_, cache.plusOpsCache, BinaryNumericOp<Op>()), false);
}
SINGLETON_NUMBER_OP_TEMPLATE(T) SingletonNumber<FT> operator+(const T& x) const {
if (x == 0) {
return *this;
}
return *this + Type(x);
}
SingletonNumber<FT>& operator+=(const Type& x) {
return *this = *this + x;
}
SINGLETON_NUMBER_OP_TEMPLATE(T) SingletonNumber<FT>& operator+=(const T& x) {
return *this = *this + x;
}
Type operator-(const Type& other) const {
if (isZero()) {
return -other;
}
if (other.isZero()) {
return *this;
}
struct Op {
FT operator()(const FT &lhs, const FT& rhs) {
return lhs - rhs;
}
};
return Type(getCachedBinaryOp(id_, other.id_, cache.minusOpsCache, BinaryNumericOp<Op>()), false);
}
SINGLETON_NUMBER_OP_TEMPLATE(T) SingletonNumber<FT> operator-(const T& x) const {
if (x == 0) {
return *this;
}
return *this - Type(x);
}
SingletonNumber<FT>& operator-=(const Type& x) {
return *this = *this - x;
}
SINGLETON_NUMBER_OP_TEMPLATE(T) SingletonNumber<FT>& operator-=(const T& x) {
return *this = *this - x;
}
Type operator*(const Type& other) const {
if (isZero()) {
return *this;
}
if (other.isZero()) {
return other;
}
if (isOne()) {
return other;
}
if (other.isOne()) {
return *this;
}
struct Op {
FT operator()(const FT &lhs, const FT& rhs) {
return lhs * rhs;
}
};
return Type(getCachedBinaryOp(id_, other.id_, cache.timesOpsCache, BinaryNumericOp<Op>()), false);
}
SINGLETON_NUMBER_OP_TEMPLATE(T) SingletonNumber<FT> operator*(const T& x) const {
if (x == 0) {
return Type(0);
}
return *this * Type(x);
}
SingletonNumber<FT>& operator*=(const Type& x) {
return *this = *this * x;
}
SINGLETON_NUMBER_OP_TEMPLATE(T) SingletonNumber<FT>& operator*=(const T& x) {
return *this = *this * x;
}
Type operator/(const Type& other) const {
if (isZero()) {
return *this;
}
if (other.isZero()) {
raise(SIGFPE); // Throw a floating point exception.
throw 0; // We won't reach this point.
}
if (other.isOne()) {
return *this;
}
struct Op {
FT operator()(const FT &lhs, const FT& rhs) {
return lhs / rhs;
}
};
return Type(getCachedBinaryOp(id_, other.id_, cache.divideOpsCache, BinaryNumericOp<Op>()), false);
}
SINGLETON_NUMBER_OP_TEMPLATE(T) SingletonNumber<FT> operator/(const T& x) const {
if (x == 0) {
return Type(0);
}
return *this / Type(x);
}
SingletonNumber<FT>& operator/=(const Type& x) {
return *this = *this / x;
}
SINGLETON_NUMBER_OP_TEMPLATE(T) SingletonNumber<FT>& operator/=(const T& x) {
return *this = *this / x;
}
};
#define SINGLETON_NUMBER_BINARY_NUMERIC_OP_FN(functionName, op) \
template <class T, class FT, std::enable_if_t<std::is_arithmetic<T>::value, bool> = true> \
SingletonNumber<FT> functionName(const T& lhs, const SingletonNumber<FT>& rhs) { \
return SingletonNumber<FT>(lhs) op rhs; \
}
#define SINGLETON_NUMBER_BINARY_BOOL_OP_FN(functionName, op) \
template <class T, class FT, std::enable_if_t<std::is_arithmetic<T>::value, bool> = true> \
bool functionName(const T& lhs, const SingletonNumber<FT>& rhs) { \
return SingletonNumber<FT>(lhs) op rhs; \
}
SINGLETON_NUMBER_BINARY_NUMERIC_OP_FN(operator+, +)
SINGLETON_NUMBER_BINARY_NUMERIC_OP_FN(operator-, -)
SINGLETON_NUMBER_BINARY_NUMERIC_OP_FN(operator*, *)
SINGLETON_NUMBER_BINARY_NUMERIC_OP_FN(operator/, /)
SINGLETON_NUMBER_BINARY_BOOL_OP_FN(operator<, <)
SINGLETON_NUMBER_BINARY_BOOL_OP_FN(operator>, >)
SINGLETON_NUMBER_BINARY_BOOL_OP_FN(operator<=, <=)
SINGLETON_NUMBER_BINARY_BOOL_OP_FN(operator>=, >=)
template<typename FT>
SingletonCache<FT> SingletonNumber<FT>::cache; // TODO define somewhere else?
template <class FT>
std::ostream& operator<<(std::ostream& out, const SingletonNumber<FT>& n) {
out << n.underlyingValue();
return out;
}
template <class FT>
std::istream& operator>>(std::istream& in, SingletonNumber<FT>& n) {
in >> n.underlyingValue();
return in;
}
namespace std {
template <class FT>
SingletonNumber<FT> sqrt(const SingletonNumber<FT>& x) {
return x.sqrt();
}
}
namespace CGAL {
template <class FT>
inline SingletonNumber<FT> min BOOST_PREVENT_MACRO_SUBSTITUTION(const SingletonNumber<FT>& x,const SingletonNumber<FT>& y){
return x.min(y);
}
template <class FT>
inline SingletonNumber<FT> max BOOST_PREVENT_MACRO_SUBSTITUTION(const SingletonNumber<FT>& x,const SingletonNumber<FT>& y){
return x.max(y);
}
CGAL_DEFINE_COERCION_TRAITS_FOR_SELF( SingletonNumber<CGAL::Gmpq>)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(short , SingletonNumber<CGAL::Gmpq>)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(int , SingletonNumber<CGAL::Gmpq>)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(long , SingletonNumber<CGAL::Gmpq>)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(float , SingletonNumber<CGAL::Gmpq>)
CGAL_DEFINE_COERCION_TRAITS_FROM_TO(double, SingletonNumber<CGAL::Gmpq>)
// CGAL_DEFINE_COERCION_TRAITS_FOR_SELF( SingletonNumber<CGAL::Gmpz>)
// CGAL_DEFINE_COERCION_TRAITS_FROM_TO(short , SingletonNumber<CGAL::Gmpz>)
// CGAL_DEFINE_COERCION_TRAITS_FROM_TO(int , SingletonNumber<CGAL::Gmpz>)
// CGAL_DEFINE_COERCION_TRAITS_FROM_TO(long , SingletonNumber<CGAL::Gmpz>)
// CGAL_DEFINE_COERCION_TRAITS_FROM_TO(float , SingletonNumber<CGAL::Gmpz>)
template <class FT>
class Algebraic_structure_traits<SingletonNumber<FT>>
: public Algebraic_structure_traits_base<SingletonNumber<FT>,Field_with_sqrt_tag>
{
typedef Algebraic_structure_traits<FT> Underlying_traits;
public:
typedef typename Underlying_traits::Is_exact Is_exact;
typedef typename Underlying_traits::Is_numerical_sensitive Is_numerical_sensitive;
};
template <class FT>
class Real_embeddable_traits<SingletonNumber<FT>>
: public INTERN_RET::Real_embeddable_traits_base<SingletonNumber<FT>, CGAL::Tag_true>
{
typedef Real_embeddable_traits<FT> Underlying_traits;
typedef SingletonNumber<FT> Type;
public:
typedef typename Underlying_traits::Is_real_embeddable Is_real_embeddable;
typedef typename Underlying_traits::Boolean Boolean;
typedef typename Underlying_traits::Comparison_result Comparison_result;
typedef typename Underlying_traits::Sign Sign;
typedef typename Underlying_traits::Is_zero Is_zero;
struct Is_finite: public CGAL::cpp98::unary_function<Type,Boolean>{
inline Boolean operator()(const Type &x) const {
return Underlying_traits::Is_finite()(x.underlyingValue());
}
};
// struct Abs: public CGAL::cpp98::unary_function<Type,Type>{
// inline Type operator()(const Type &x) const {
// return Underlying_traits::Abs()(x.underlyingValue());
// }
// };
struct To_interval: public CGAL::cpp98::unary_function<Type,std::pair<double,double> >{
inline std::pair<double,double> operator()(const Type &x) const {
typename Underlying_traits::To_interval to_interval;
return to_interval(x.underlyingValue());
}
};
};
template <>
class Fraction_traits< SingletonNumber<CGAL::Gmpq> > {
public:
typedef SingletonNumber<CGAL::Gmpq> Type;
typedef ::CGAL::Tag_true Is_fraction;
typedef Gmpz Numerator_type;
typedef Gmpz Denominator_type;
typedef Algebraic_structure_traits< Gmpz >::Gcd Common_factor;
class Decompose {
public:
typedef SingletonNumber<CGAL::Gmpq> first_argument_type;
typedef Gmpz& second_argument_type;
typedef Gmpz& third_argument_type;
void operator () (const SingletonNumber<CGAL::Gmpq>& rat, Gmpz& num,Gmpz& den) {
// TODO(ochafik): compose proper singleton num and denom, and their transition to a singleton rational.
num = rat.underlyingValue().numerator();
den = rat.underlyingValue().denominator();
}
};
class Compose {
public:
typedef Gmpz first_argument_type;
typedef Gmpz second_argument_type;
typedef SingletonNumber<CGAL::Gmpq> result_type;
SingletonNumber<CGAL::Gmpq> operator () (const Gmpz& num,const Gmpz& den) {
// TODO(ochafik): compose proper singleton num and denom, and their transition to a singleton rational.
return SingletonNumber<CGAL::Gmpq>(CGAL::Gmpq(num, den));
}
};
};
class SingletonNumberGMP_arithmetic_kernel : public internal::Arithmetic_kernel_base {
public:
typedef Gmpz Integer;
typedef SingletonNumber<CGAL::Gmpq> Rational;
};
template <>
struct Get_arithmetic_kernel<SingletonNumber<Gmpz>> {
typedef SingletonNumberGMP_arithmetic_kernel Arithmetic_kernel;
};
template <>
struct Get_arithmetic_kernel<SingletonNumber<Gmpq>> {
typedef SingletonNumberGMP_arithmetic_kernel Arithmetic_kernel;
};
} //namespace CGAL
namespace Eigen {
template<class> struct NumTraits;
template<> struct NumTraits<SingletonNumber<CGAL::Gmpq>>
{
typedef SingletonNumber<CGAL::Gmpq> Real;
typedef SingletonNumber<CGAL::Gmpq> NonInteger;
typedef SingletonNumber<CGAL::Gmpq> Nested;
typedef SingletonNumber<CGAL::Gmpq> Literal;
static inline Real epsilon() { return 0; }
static inline Real dummy_precision() { return 0; }
enum {
IsInteger = 0,
IsSigned = 1,
IsComplex = 0,
RequireInitialization = 1,
ReadCost = 6,
AddCost = 150, // TODO(ochafik): reduce this
MulCost = 100 // TODO(ochafik): reduce this
};
};
namespace internal {
template<>
struct significant_decimals_impl<SingletonNumber<CGAL::Gmpq>>
{
static inline int run()
{
return 0;
}
};
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment