0

Consider the following minimal example:

#include <stdlib.h>

#include "CGAL\Exact_predicates_exact_constructions_kernel.h"

typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::FT NumberType;

int main()
{
    int numberOfIterations = 4000;

    NumberType sum = 0;
    for (int i = 0; i < numberOfIterations; i++)
    {
        sum += 1;
    }

    std::cout << "Sum = " << sum << std::endl;

    getchar();

    return 0;
}

All seems to work fine if the numberOfIterations parameter assumes small values. However, when it is set to a value greater than some thousands, the program outputs the correct value for sum, but gives a Stack overflow error (parameters: 0x0000000000000001, 0x000000C4D2EF3FF8) when deleting the memory (I suppose).

I read the CGAL documentation on Number Types, but it is not clear to me the reason why this happens. What should I do to perform operations on exact numbers in a correct way?

Here is a more complicated example (as requested in the comments below), strictly related to my original problem:

#include <stdlib.h>

#include <CGAL/number_utils.h>
#include "CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h"
#include "CGAL/Polyhedron_3.h"

typedef CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt Kernel;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3;
typedef Kernel::FT NumberType;
typedef Kernel::Vector_3 Vector_3;
typedef Kernel::Direction_3 Direction_3;

Direction_3 sampleFunction(std::vector<Vector_3> vectors)
{
    NumberType denominator = 0;
    NumberType numeratorX = 0;
    NumberType numeratorY = 0;
    NumberType numeratorZ = 0;

    for (int vf = 0; vf < vectors.size(); vf++)
    {
        Vector_3 vectorF = vectors.at(vf);

        numeratorX += vectorF.x();
        numeratorY += vectorF.y();
        numeratorZ += vectorF.z();

        for (int vg = 0; vg < vectors.size(); vg++)
        {
            Vector_3 vectorG = vectors.at(vg);

            denominator += CGAL::scalar_product(vectorF, vectorG);
        }
    }
    NumberType lambda = CGAL::sqrt(denominator);

    NumberType x = numeratorX / lambda;
    NumberType y = numeratorY / lambda;
    NumberType z = numeratorZ / lambda;

    Vector_3 result(x, y, z);

    return result.direction();
}

int main()
{
    int numberOfVectors = 100;

    std::vector<Vector_3> vectors;
    srand(0);
    for (int i = 0; i < numberOfVectors; i++)
    {
        Vector_3 newVector(rand() % 100 - 50, rand() % 100 - 50, rand() % 100 - 50);
        newVector /= CGAL::sqrt(newVector.squared_length());
        vectors.push_back(newVector);
    }

    Direction_3 direction = sampleFunction(vectors);

    std::cout << direction.dx() << " " << direction.dy() << " " << direction.dz() << " " << std::endl;

    return 0;
}
Vadim Kotov
  • 8,084
  • 8
  • 48
  • 62
simonet
  • 295
  • 1
  • 14
  • With your configuration, FT is probably Lazy_exact_nt, which builds a tree that remembers all the operations performed on it. Sometimes evaluation has to go through this tree using a recursive function, which requires a large stack. So does destruction (though I thought I saw either an issue or a pull request about that). – Marc Glisse Nov 22 '18 at 11:47
  • https://github.com/CGAL/cgal/issues/1118 . Something like this branch could also help: https://github.com/mglisse/cgal/tree/Number_types-lazy-glisse . – Marc Glisse Nov 22 '18 at 11:48
  • Thanks for the explanation and the links. Let's see if I understood: (1) the value of `sum` that is printed is only an approximated value (otherwise the computation would lead to a stack overflow already there). (2) there is no way (at the moment) for doing a large number of operations on exact numbers (except by downloading/installing/linking to the branch you have provided). Am I right? – simonet Nov 23 '18 at 13:41
  • (1) yes but as long as it fits in double the approximation is perfect. (2) you can call sum.exact() once in a while to collapse the tree, or better use a non-lazy type. You can also increase the size of the stack (there are many questions on this site explaining how on various systems). – Marc Glisse Nov 23 '18 at 15:03
  • Your example code is too simple to be useful. There is a category of computations, iterated/cascaded computations, for which this CGAL number type is not appropriate in general, even though it can be tweaked to make it work for your particular simple case, as Marc Glisse mentions. Do you have a more "real world" example of thing you would like to do, for which this is failing as well ? – Sylvain Pion Nov 23 '18 at 15:29
  • @MarcGlisse Thanks for the clarifications. So one option is to use a different exact kernel. I tried `Exact_predicates_exact_constructions_kernel_with_sqrt`, that uses `CORE::Expr` as number type, but it resulted in a stack overflow occurred during the `std::cout` command. – simonet Nov 23 '18 at 17:26
  • @SylvainPion I provided a more specific example at the bottom of the post. – simonet Nov 23 '18 at 17:27
  • I had forgotten that CORE::Expr was building trees as well... You could go with `CGAL::Simple_cartesian` for something simple. – Marc Glisse Nov 23 '18 at 18:38
  • I didn't read your last example carefully, but do you really need exact constructions for it? – Marc Glisse Nov 23 '18 at 18:40
  • Ok, yes... If I wanted to solve the problem (`sampleFunction` computes the result of a formula obtained with the _method of Lagrange Multipliers_), I should use an exact kernel. In the real application, though, I think I do not really need exact precision. – simonet Nov 26 '18 at 07:25

0 Answers0