diff --git a/ChangeLog.txt b/ChangeLog.txt index d4018b7af..2679aa0dd 100644 --- a/ChangeLog.txt +++ b/ChangeLog.txt @@ -20,6 +20,7 @@ - Fix segmenting a simple curve. - Fix export tiled pdf on Mac OS. - Improve visualization for tools. +- Fix calculating an elliptical arc. # Valentina 0.7.51 April 18, 2022 - Z value change for a layout piece. diff --git a/src/app/puzzle/carousel/vpcarrousel.h b/src/app/puzzle/carousel/vpcarrousel.h index 3f160cafc..ed3ad74b4 100644 --- a/src/app/puzzle/carousel/vpcarrousel.h +++ b/src/app/puzzle/carousel/vpcarrousel.h @@ -77,7 +77,6 @@ public: * @brief Clear Clears the carrousel (removes everything) */ void Clear(); - auto Layout() const -> VPLayoutWeakPtr; public slots: diff --git a/src/libs/vgeometry/vellipticalarc.cpp b/src/libs/vgeometry/vellipticalarc.cpp index 0811168da..7554aae4f 100644 --- a/src/libs/vgeometry/vellipticalarc.cpp +++ b/src/libs/vgeometry/vellipticalarc.cpp @@ -37,9 +37,171 @@ #include "../ifc/ifcdef.h" #include "../ifc/exception/vexception.h" #include "../vmisc/vabstractapplication.h" +#include "../vmisc/fpm/fixed.hpp" +#include "../vmisc/fpm/math.hpp" #include "../vmisc/compatibility.h" #include "vabstractcurve.h" #include "vellipticalarc_p.h" +#include "../vmisc/vmath.h" + +namespace +{ +constexpr qreal tolerance = accuracyPointOnLine/32; +//--------------------------------------------------------------------------------------------------------------------- +auto VLen(fpm::fixed_16_16 x, fpm::fixed_16_16 y) -> fpm::fixed_16_16 +{ + x = fpm::abs(x); + y = fpm::abs(y); + if (x > y) + { + return x + qMax(y/8, y/2 - x/8); + } + + return y + qMax(x/8, x/2 - y/8); +} + +//--------------------------------------------------------------------------------------------------------------------- +auto AuxRadius(fpm::fixed_16_16 xP, fpm::fixed_16_16 yP, fpm::fixed_16_16 xQ, fpm::fixed_16_16 yQ) -> fpm::fixed_16_16 +{ + fpm::fixed_16_16 dP = VLen(xP, yP); + fpm::fixed_16_16 dQ = VLen(xQ, yQ); + fpm::fixed_16_16 dJ = VLen(xP + xQ, yP + yQ); + fpm::fixed_16_16 dK = VLen(xP - xQ, yP - yQ); + fpm::fixed_16_16 r1 = qMax(dP, dQ); + fpm::fixed_16_16 r2 = qMax(dJ, dK); + return qMax(r1 + r1/16, r2 - r2/4); +} + +//--------------------------------------------------------------------------------------------------------------------- +auto AngularInc(fpm::fixed_16_16 xP, fpm::fixed_16_16 yP, fpm::fixed_16_16 xQ, fpm::fixed_16_16 yQ, + fpm::fixed_16_16 flatness) -> int +{ + + fpm::fixed_16_16 r = AuxRadius(xP, yP, xQ, yQ); + fpm::fixed_16_16 err2{r >> 3}; + // 2nd-order term + fpm::fixed_16_16 err4{r >> 7}; + // 4th-order term + const int kmax = qRound(0.5 * std::log2(maxSceneSize / (8. * tolerance))); + for (int k = 0; k < kmax; ++k) + { + if (flatness >= err2 + err4) + { + return k; + } + err2 >>= 2; + err4 >>= 4; + } + return kmax; +} + +//--------------------------------------------------------------------------------------------------------------------- +inline void CircleGen(fpm::fixed_16_16& u, fpm::fixed_16_16& v, uint k) +{ + u -= v >> k; + v += u >> k; +} + +//--------------------------------------------------------------------------------------------------------------------- +auto InitialValue(fpm::fixed_16_16 u0, fpm::fixed_16_16 v0, uint k) -> fpm::fixed_16_16 +{ + uint shift = 2*k + 3; + fpm::fixed_16_16 w {u0 >> shift}; + + fpm::fixed_16_16 U0 = u0 - w + (v0 >> (k + 1)); + w >>= (shift + 1); + U0 -= w; + w >>= shift; + U0 -= w; + return U0; +} + +//--------------------------------------------------------------------------------------------------------------------- +auto EllipseCore(fpm::fixed_16_16 xC, fpm::fixed_16_16 yC, fpm::fixed_16_16 xP, fpm::fixed_16_16 yP, + fpm::fixed_16_16 xQ, fpm::fixed_16_16 yQ, fpm::fixed_16_16 sweep, + fpm::fixed_16_16 flatness) -> QVector +{ + uint k = qMin(static_cast(AngularInc(xP, yP, xQ, yQ, flatness)), 16U); + const uint count = static_cast(sweep.raw_value()) >> (16 - k); + + QVector arc; + arc.reserve(static_cast(count) + 1); + + // Arc start point + arc.append({static_cast(xP + xC), static_cast(yP + yC)}); + + xQ = InitialValue(xQ, xP, k); + yQ = InitialValue(yQ, yP, k); + + for (uint i = 0; i < count; ++i) + { + CircleGen(xQ, xP, k); + CircleGen(yQ, yP, k); + arc.append({static_cast(xP + xC), static_cast(yP + yC)}); + } + + return arc; +} + +//--------------------------------------------------------------------------------------------------------------------- +auto EllipticArcPoints(QPointF c, qreal radius1, qreal radius2, qreal astart, qreal asweep, + qreal approximationScale) -> QVector +{ + fpm::fixed_16_16 xC{c.x()}; + fpm::fixed_16_16 yC{c.y()}; + + fpm::fixed_16_16 xP{c.x() + radius1}; + fpm::fixed_16_16 yP{c.y()}; + + fpm::fixed_16_16 xQ{c.x()}; + fpm::fixed_16_16 yQ{c.y() - radius2}; + + xP -= xC; + yP -= yC; + xQ -= xC; + yQ -= yC; + + if (not qFuzzyIsNull(astart)) + { + // Set new conjugate diameter end points P’ and Q’ + fpm::fixed_16_16 cosa {cos(astart)}; + fpm::fixed_16_16 sina {sin(astart)}; + fpm::fixed_16_16 x {xP * cosa + xQ * sina}; + fpm::fixed_16_16 y {yP * cosa + yQ * sina}; + + xQ = xQ * cosa - xP * sina; + yQ = yQ * cosa - yP * sina; + xP = x; + yP = y; + } + + // If sweep angle is negative, switch direction + if (asweep < 0) + { + xQ = -xQ; + yQ = -yQ; + asweep = -asweep; + } + + if(approximationScale < minCurveApproximationScale || approximationScale > maxCurveApproximationScale) + { + approximationScale = VAbstractApplication::VApp()->Settings()->GetCurveApproximationScale(); + } + + fpm::fixed_16_16 flatness {maxCurveApproximationScale / approximationScale * tolerance}; + fpm::fixed_16_16 swangle{asweep}; + QVector arc = EllipseCore(xC, yC, xP, yP, xQ, yQ, swangle, flatness); + + // Arc end point + fpm::fixed_16_16 cosb {qCos(asweep)}; + fpm::fixed_16_16 sinb {qSin(asweep)}; + xP = xP*cosb + xQ*sinb; + yP = yP*cosb + yQ*sinb; + arc.append({static_cast(xP+xC), static_cast(yP+yC)}); + + return arc; +} +} // namespace //--------------------------------------------------------------------------------------------------------------------- /** @@ -59,11 +221,11 @@ VEllipticalArc::VEllipticalArc() * @param f2 end angle (degree). */ VEllipticalArc::VEllipticalArc (const VPointF ¢er, qreal radius1, qreal radius2, const QString &formulaRadius1, - const QString &formulaRadius2, qreal f1, const QString &formulaF1, qreal f2, - const QString &formulaF2, qreal rotationAngle, const QString &formulaRotationAngle, - quint32 idObject, Draw mode) + const QString &formulaRadius2, qreal f1, const QString &formulaF1, qreal f2, + const QString &formulaF2, qreal rotationAngle, const QString &formulaRotationAngle, + quint32 idObject, Draw mode) : VAbstractArc(GOType::EllipticalArc, center, f1, formulaF1, f2, formulaF2, idObject, mode), - d (new VEllipticalArcData(radius1, radius2, formulaRadius1, formulaRadius2, rotationAngle, formulaRotationAngle)) + d (new VEllipticalArcData(radius1, radius2, formulaRadius1, formulaRadius2, rotationAngle, formulaRotationAngle)) { CreateName(); } @@ -72,7 +234,7 @@ VEllipticalArc::VEllipticalArc (const VPointF ¢er, qreal radius1, qreal radi VEllipticalArc::VEllipticalArc(const VPointF ¢er, qreal radius1, qreal radius2, qreal f1, qreal f2, qreal rotationAngle) : VAbstractArc(GOType::EllipticalArc, center, f1, f2, NULL_ID, Draw::Calculation), - d (new VEllipticalArcData(radius1, radius2, rotationAngle)) + d (new VEllipticalArcData(radius1, radius2, rotationAngle)) { CreateName(); } @@ -83,7 +245,7 @@ VEllipticalArc::VEllipticalArc(qreal length, const QString &formulaLength, const const QString &formulaF1, qreal rotationAngle, const QString &formulaRotationAngle, quint32 idObject, Draw mode) : VAbstractArc(GOType::EllipticalArc, formulaLength, center, f1, formulaF1, idObject, mode), - d (new VEllipticalArcData(radius1, radius2, formulaRadius1, formulaRadius2, rotationAngle, formulaRotationAngle)) + d (new VEllipticalArcData(radius1, radius2, formulaRadius1, formulaRadius2, rotationAngle, formulaRotationAngle)) { CreateName(); FindF2(length); @@ -93,7 +255,7 @@ VEllipticalArc::VEllipticalArc(qreal length, const QString &formulaLength, const VEllipticalArc::VEllipticalArc(qreal length, const VPointF ¢er, qreal radius1, qreal radius2, qreal f1, qreal rotationAngle) : VAbstractArc(GOType::EllipticalArc, center, f1, NULL_ID, Draw::Calculation), - d (new VEllipticalArcData(radius1, radius2, rotationAngle)) + d (new VEllipticalArcData(radius1, radius2, rotationAngle)) { CreateName(); FindF2(length); @@ -193,7 +355,7 @@ auto VEllipticalArc::Move(qreal length, qreal angle, const QString &prefix) cons const VPointF center = oldCenter.Move(length, angle); const QPointF position = d->m_transform.inverted().map(center.toQPointF()) - - d->m_transform.inverted().map(oldCenter.toQPointF()); + d->m_transform.inverted().map(oldCenter.toQPointF()); QTransform t = d->m_transform; t.translate(position.x(), position.y()); @@ -285,45 +447,19 @@ auto VEllipticalArc::GetCenter() const -> VPointF auto VEllipticalArc::GetPoints() const -> QVector { const QPointF center = VAbstractArc::GetCenter().toQPointF(); - QRectF box(center.x() - d->radius1, center.y() - d->radius2, d->radius1*2, d->radius2*2); - QLineF startLine(center.x(), center.y(), center.x() + d->radius1, center.y()); - QLineF endLine = startLine; - - startLine.setAngle(VAbstractArc::GetStartAngle()); - endLine.setAngle(RealEndAngle()); - qreal sweepAngle = startLine.angleTo(endLine); - - if (qFuzzyIsNull(sweepAngle)) - { - sweepAngle = 360; - } - - QPainterPath path; - path.moveTo(GetP1()); - path.arcTo(box, VAbstractArc::GetStartAngle(), sweepAngle); - path.moveTo(GetP2()); + // Generate complete ellipse because angles are not correct and have to be fixed manually + QVector points = EllipticArcPoints(center, d->radius1, d->radius2, 0.0, M_2PI, GetApproximationScale()); + points = ArcPoints(points); QTransform t = d->m_transform; t.translate(center.x(), center.y()); t.rotate(-GetRotationAngle()); t.translate(-center.x(), -center.y()); - path = t.map(path); + std::transform(points.begin(), points.end(), points.begin(), [t](const QPointF &point) { return t.map(point); }); - QPolygonF polygon; - const QList sub = path.toSubpathPolygons(); - if (not sub.isEmpty()) - { - polygon = ConstFirst (path.toSubpathPolygons()); - - if (not polygon.isEmpty() && not VFuzzyComparePoints(GetP1(), ConstFirst (polygon))) - { - polygon.removeFirst(); // remove point (0;0) - } - } - - return static_cast>(polygon); + return points; } //--------------------------------------------------------------------------------------------------------------------- @@ -360,7 +496,7 @@ auto VEllipticalArc::CutArc(const qreal &length, VEllipticalArc &arc1, VElliptic const QString errorMsg = QObject::tr("Unable to cut curve '%1'. The curve is too short.").arg(name()); VAbstractApplication::VApp()->IsPedantic() ? throw VException(errorMsg) : - qWarning() << VAbstractApplication::warningMessageSignature + errorMsg; + qWarning() << VAbstractApplication::warningMessageSignature + errorMsg; return {}; } @@ -380,10 +516,10 @@ auto VEllipticalArc::CutArc(const qreal &length, VEllipticalArc &arc1, VElliptic else { errorMsg = QObject::tr("Curve '%1'. Length of a cut segment is too small. Optimize it to minimal value.") - .arg(name()); + .arg(name()); } VAbstractApplication::VApp()->IsPedantic() ? throw VException(errorMsg) : - qWarning() << VAbstractApplication::warningMessageSignature + errorMsg; + qWarning() << VAbstractApplication::warningMessageSignature + errorMsg; } else if (length > maxLength) { @@ -393,15 +529,15 @@ auto VEllipticalArc::CutArc(const qreal &length, VEllipticalArc &arc1, VElliptic if (not pointName.isEmpty()) { errorMsg = QObject::tr("Curve '%1'. Length of a cut segment (%2) is too big. Optimize it to maximal value.") - .arg(name(), pointName); + .arg(name(), pointName); } else { errorMsg = QObject::tr("Curve '%1'. Length of a cut segment is too big. Optimize it to maximal value.") - .arg(name()); + .arg(name()); } VAbstractApplication::VApp()->IsPedantic() ? throw VException(errorMsg) : - qWarning() << VAbstractApplication::warningMessageSignature + errorMsg; + qWarning() << VAbstractApplication::warningMessageSignature + errorMsg; } else { @@ -410,13 +546,13 @@ auto VEllipticalArc::CutArc(const qreal &length, VEllipticalArc &arc1, VElliptic // the first arc has given length and startAngle just like in the origin arc arc1 = VEllipticalArc (len, QString().setNum(length), GetCenter(), d->radius1, d->radius2, - d->formulaRadius1, d->formulaRadius2, GetStartAngle(), GetFormulaF1(), d->rotationAngle, - GetFormulaRotationAngle(), getIdObject(), getMode()); + d->formulaRadius1, d->formulaRadius2, GetStartAngle(), GetFormulaF1(), d->rotationAngle, + GetFormulaRotationAngle(), getIdObject(), getMode()); // the second arc has startAngle just like endAngle of the first arc // and it has endAngle just like endAngle of the origin arc arc2 = VEllipticalArc (GetCenter(), d->radius1, d->radius2, d->formulaRadius1, d->formulaRadius2, - arc1.GetEndAngle(), arc1.GetFormulaF2(), GetEndAngle(), GetFormulaF2(), d->rotationAngle, - GetFormulaRotationAngle(), getIdObject(), getMode()); + arc1.GetEndAngle(), arc1.GetFormulaF2(), GetEndAngle(), GetFormulaF2(), d->rotationAngle, + GetFormulaRotationAngle(), getIdObject(), getMode()); return arc1.GetP1(); } @@ -432,7 +568,7 @@ auto VEllipticalArc::CutArc(const qreal &length, const QString &pointName) const //--------------------------------------------------------------------------------------------------------------------- void VEllipticalArc::CreateName() { - QString name = ELARC_ + QString("%1").arg(this->GetCenter().name()); + QString name = ELARC_ + QStringLiteral("%1").arg(this->GetCenter().name()); const QString nameStr = QStringLiteral("_%1"); if (getMode() == Draw::Modeling && getIdObject() != NULL_ID) @@ -555,23 +691,101 @@ auto VEllipticalArc::GetP(qreal angle) const -> QPointF } //--------------------------------------------------------------------------------------------------------------------- -auto VEllipticalArc::RealEndAngle() const -> qreal +auto VEllipticalArc::ArcPoints(QVector points) const -> QVector { - qreal endAngle = VEllipticalArc::OptimizeAngle(VAbstractArc::GetEndAngle()); - - if (qFuzzyIsNull(endAngle) || - VFuzzyComparePossibleNulls(endAngle, 90) || - VFuzzyComparePossibleNulls(endAngle, 180) || - VFuzzyComparePossibleNulls(endAngle, 270) || - VFuzzyComparePossibleNulls(endAngle, 360)) + if (points.size() < 2 || VFuzzyComparePossibleNulls(VAbstractArc::GetStartAngle(), VAbstractArc::GetEndAngle())) { - return endAngle; + return points; } - endAngle = qRadiansToDegrees(qAtan2(d->radius1 * qSin(qDegreesToRadians(endAngle)), - d->radius2 * qCos(qDegreesToRadians(endAngle)))); + QPointF center = static_cast(GetCenter()); + qreal radius = qMax(d->radius1, d->radius2) * 2; - return endAngle; + QLineF start(center.x(), center.y(), center.x() + radius, center.y()); + start.setAngle(VAbstractArc::GetStartAngle()); + + QLineF end(center.x(), center.y(), center.x() + radius, center.y()); + end.setAngle(VAbstractArc::GetEndAngle()); + + if (VFuzzyComparePossibleNulls(start.angle(), end.angle())) + { + return points; + } + + auto IsBoundedIntersection = [](QLineF::IntersectType type, QPointF p, const QLineF &segment1, + const QLineF &segment2) + { + return type == QLineF::BoundedIntersection || + (VGObject::IsPointOnLineSegment (p, segment1.p1(), segment2.p1()) && + VGObject::IsPointOnLineSegment (p, segment2.p1(), segment2.p2())); + }; + + bool begin = true; + + if (start.angle() > end.angle()) + { + for (int i=0; i < points.size()-1; ++i) + { + QLineF edge(points.at(i), points.at(i+1)); + + QPointF p; + QLineF::IntersectType type = Intersects(start, edge, &p); + + // QLineF::intersects not always accurate on edge cases + if (IsBoundedIntersection(type, p, edge, start)) + { + QVector head = points.mid(0, i+1); + QVector tail = points.mid(i+1, -1); + tail.prepend(p); + points = tail + head; + begin = false; + break; + } + } + } + + QVector arc; + arc.reserve(points.size()); + + for (int i=0; i < points.size()-1; ++i) + { + QLineF edge(points.at(i), points.at(i+1)); + + if (begin) + { + QPointF p; + QLineF::IntersectType type = Intersects(start, edge, &p); + + // QLineF::intersects not always accurate on edge cases + if (IsBoundedIntersection(type, p, edge, start)) + { + arc.append(p); + begin = false; + } + } + else + { + QPointF p; + QLineF::IntersectType type = Intersects(end, edge, &p); + + // QLineF::intersects not always accurate on edge cases + if (IsBoundedIntersection(type, p, edge, end)) + { + arc.append(points.at(i)); + arc.append(p); + break; + } + + arc.append(points.at(i)); + } + } + + if (arc.isEmpty()) + { + return points; + } + + return arc; } //--------------------------------------------------------------------------------------------------------------------- diff --git a/src/libs/vgeometry/vellipticalarc.h b/src/libs/vgeometry/vellipticalarc.h index 11ba1e17e..c69148e61 100644 --- a/src/libs/vgeometry/vellipticalarc.h +++ b/src/libs/vgeometry/vellipticalarc.h @@ -114,10 +114,8 @@ private: QSharedDataPointer d; auto MaxLength() const -> qreal; - auto GetP(qreal angle) const -> QPointF; - - auto RealEndAngle() const -> qreal; + auto ArcPoints(QVector points) const -> QVector; }; Q_DECLARE_METATYPE(VEllipticalArc) // NOLINT diff --git a/src/libs/vmisc/def.cpp b/src/libs/vmisc/def.cpp index 361fc7f33..286937611 100644 --- a/src/libs/vmisc/def.cpp +++ b/src/libs/vmisc/def.cpp @@ -64,13 +64,6 @@ #include "../ifc/exception/vexception.h" #include "literals.h" -const qreal defCurveApproximationScale = 0.5; -const qreal minCurveApproximationScale = 0.2; -const qreal maxCurveApproximationScale = 10.0; - -const int minLabelFontSize = 5; -const int maxLabelFontSize = 100; - //--------------------------------------------------------------------------------------------------------------------- QPixmap QPixmapFromCache(const QString &pixmapPath) { diff --git a/src/libs/vmisc/def.h b/src/libs/vmisc/def.h index 353622c0f..3d2d3cb8e 100644 --- a/src/libs/vmisc/def.h +++ b/src/libs/vmisc/def.h @@ -74,13 +74,14 @@ class QMarginsF; class VTranslateMeasurements; class QGraphicsItem; -#define SceneSize 50000 -extern const qreal defCurveApproximationScale; -extern const qreal minCurveApproximationScale; -extern const qreal maxCurveApproximationScale; +constexpr qreal maxSceneSize = ((20.0 * 1000.0) / 25.4) * PrintDPI; // 20 meters in pixels -extern const int minLabelFontSize; -extern const int maxLabelFontSize; +constexpr qreal defCurveApproximationScale = 0.5; +constexpr qreal minCurveApproximationScale = 0.2; +constexpr qreal maxCurveApproximationScale = 10.0; + +constexpr int minLabelFontSize = 5; +constexpr int maxLabelFontSize = 100; enum class NodeDetail : qint8 { Contour, Modeling }; enum class SceneObject : qint8 { Point, Line, Spline, Arc, ElArc, SplinePath, Detail, Unknown }; diff --git a/src/libs/vmisc/fpm/LICENSE b/src/libs/vmisc/fpm/LICENSE new file mode 100644 index 000000000..bb86b71a8 --- /dev/null +++ b/src/libs/vmisc/fpm/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2019 Mike Lankamp + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/src/libs/vmisc/fpm/fixed.hpp b/src/libs/vmisc/fpm/fixed.hpp new file mode 100644 index 000000000..5471a0dc3 --- /dev/null +++ b/src/libs/vmisc/fpm/fixed.hpp @@ -0,0 +1,497 @@ +#ifndef FPM_FIXED_HPP +#define FPM_FIXED_HPP + +#include +#include +#include +#include +#include +#include +#include + +#if QT_VERSION < QT_VERSION_CHECK(5, 5, 0) +#include "../diagnostic.h" +#endif // QT_VERSION < QT_VERSION_CHECK(5, 5, 0) + +namespace fpm +{ + +//! Fixed-point number type +//! \tparam BaseType the base integer type used to store the fixed-point number. This can be a signed or +//! unsigned type. +//! \tparam IntermediateType the integer type used to store intermediate results during calculations. +//! \tparam FractionBits the number of bits of the BaseType used to store the fraction +template +class fixed +{ + static_assert(std::is_integral::value, "BaseType must be an integral type"); + static_assert(FractionBits > 0, "FractionBits must be greater than zero"); + static_assert(FractionBits <= sizeof(BaseType) * 8 - 1, "BaseType must at least be able to contain entire " + "fraction, with space for at least one integral bit"); + static_assert(sizeof(IntermediateType) > sizeof(BaseType), "IntermediateType must be larger than BaseType"); + static_assert(std::is_signed::value == std::is_signed::value, + "IntermediateType must have same signedness as BaseType"); + + // Although this value fits in the BaseType in terms of bits, if there's only one integral bit, this value + // is incorrect (flips from positive to negative), so we must extend the size to IntermediateType. + static constexpr IntermediateType FRACTION_MULT = IntermediateType(1) << FractionBits; + + struct raw_construct_tag {}; + constexpr inline fixed(BaseType val, raw_construct_tag /*unused*/) noexcept : m_value(val) {} + +public: + inline fixed() noexcept = default; + + // Converts an integral number to the fixed-point type. + // Like static_cast, this truncates bits that don't fit. + template ::value>::type* = nullptr> + constexpr inline explicit fixed(T val) noexcept + : m_value(static_cast(val * FRACTION_MULT)) + {} + + // Converts an floating-point number to the fixed-point type. + // Like static_cast, this truncates bits that don't fit. + template ::value>::type* = nullptr> + constexpr inline explicit fixed(T val) noexcept + : m_value(static_cast((val >= 0.0) ? (val * FRACTION_MULT + T{0.5}) : (val * FRACTION_MULT - T{0.5}))) + {} + + // Constructs from another fixed-point type with possibly different underlying representation. + // Like static_cast, this truncates bits that don't fit. + template + constexpr inline explicit fixed(fixed val) noexcept + : m_value(from_fixed_point(val.raw_value()).raw_value()) + {} + + // Explicit conversion to a floating-point type + template ::value>::type* = nullptr> + constexpr inline explicit operator T() const noexcept + { + return static_cast(m_value) / FRACTION_MULT; + } + + // Explicit conversion to an integral type + template ::value>::type* = nullptr> + constexpr inline explicit operator T() const noexcept + { + return static_cast(m_value / FRACTION_MULT); + } + + // Returns the raw underlying value of this type. + // Do not use this unless you know what you're doing. + constexpr inline auto raw_value() const noexcept -> BaseType + { + return m_value; + } + + //! Constructs a fixed-point number from another fixed-point number. + //! \tparam NumFractionBits the number of bits used by the fraction in \a value. + //! \param value the integer fixed-point number + template FractionBits)>::type* = nullptr> + static constexpr inline auto from_fixed_point(T value) noexcept -> fixed + { + // To correctly round the last bit in the result, we need one more bit of information. + // We do this by multiplying by two before dividing and adding the LSB to the real result. + return fixed(static_cast( + value / (T(1) << (NumFractionBits - FractionBits)) + + (value / (T(1) << (NumFractionBits - FractionBits - 1)) % 2)), + raw_construct_tag{}); + } + + template ::type* = nullptr> + static constexpr inline auto from_fixed_point(T value) noexcept -> fixed + { + return fixed(static_cast( + value * (T(1) << (FractionBits - NumFractionBits))), + raw_construct_tag{}); + } + + // Constructs a fixed-point number from its raw underlying value. + // Do not use this unless you know what you're doing. + static constexpr inline auto from_raw_value(BaseType value) noexcept -> fixed + { + return fixed(value, raw_construct_tag{}); + } + + // + // Constants + // + static constexpr auto e() -> fixed { return from_fixed_point<61>(6267931151224907085LL); } + static constexpr auto pi() -> fixed { return from_fixed_point<61>(7244019458077122842LL); } + static constexpr auto half_pi() -> fixed { return from_fixed_point<62>(7244019458077122842LL); } + static constexpr auto two_pi() -> fixed { return from_fixed_point<60>(7244019458077122842LL); } + + // + // Arithmetic member operators + // + + constexpr inline auto operator-() const noexcept -> fixed + { + return fixed::from_raw_value(-m_value); + } + + inline auto operator+=(const fixed& y) noexcept -> fixed& + { + m_value += y.m_value; + return *this; + } + + template ::value>::type* = nullptr> + inline auto operator+=(I y) noexcept -> fixed& + { + m_value += y * FRACTION_MULT; + return *this; + } + + inline auto operator-=(const fixed& y) noexcept -> fixed& + { + m_value -= y.m_value; + return *this; + } + + template ::value>::type* = nullptr> + inline auto operator-=(I y) noexcept -> fixed& + { + m_value -= y * FRACTION_MULT; + return *this; + } + + inline auto operator*=(const fixed& y) noexcept -> fixed& + { + // Normal fixed-point multiplication is: x * y / 2**FractionBits. + // To correctly round the last bit in the result, we need one more bit of information. + // We do this by multiplying by two before dividing and adding the LSB to the real result. + auto value = static_cast( + static_cast(static_cast(m_value) * y.m_value) / (FRACTION_MULT / 2)); + m_value = static_cast((value / 2) + (value % 2)); + return *this; + } + + template ::value>::type* = nullptr> + inline auto operator*=(I y) noexcept -> fixed& + { + m_value *= y; + return *this; + } + + inline auto operator/=(const fixed& y) noexcept -> fixed& + { + assert(y.m_value != 0); + // Normal fixed-point division is: x * 2**FractionBits / y. + // To correctly round the last bit in the result, we need one more bit of information. + // We do this by multiplying by two before dividing and adding the LSB to the real result. + auto value = (static_cast(m_value) * FRACTION_MULT * 2) / y.m_value; + m_value = static_cast((value / 2) + (value % 2)); + return *this; + } + + template ::value>::type* = nullptr> + inline auto operator/=(I y) noexcept -> fixed& + { + m_value /= y; + return *this; + } + + // + // Shift operators + // + + template ::value>::type* = nullptr> + inline auto operator>>=(I y) noexcept -> fixed& + { + m_value >>= y; + return *this; + } + + template ::value>::type* = nullptr> + inline auto operator<<=(I y) noexcept -> fixed& + { + m_value <<= y; + return *this; + } + +private: + BaseType m_value; +}; + +// +// Convenience typedefs +// + +using fixed_16_16 = fixed; +using fixed_24_8 = fixed; +using fixed_8_24 = fixed; + +// +// Addition +// + +template +constexpr inline auto operator+(const fixed& x, const fixed& y) noexcept -> fixed +{ + return fixed(x) += y; +} + +template ::value>::type* = nullptr> +constexpr inline auto operator+(const fixed& x, T y) noexcept -> fixed +{ + return fixed(x) += y; +} + +template ::value>::type* = nullptr> +constexpr inline auto operator+(T x, const fixed& y) noexcept -> fixed +{ + return fixed(y) += x; +} + +// +// Subtraction +// + +template +constexpr inline auto operator-(const fixed& x, const fixed& y) noexcept -> fixed +{ + return fixed(x) -= y; +} + +template ::value>::type* = nullptr> +constexpr inline auto operator-(const fixed& x, T y) noexcept -> fixed +{ + return fixed(x) -= y; +} + +template ::value>::type* = nullptr> +constexpr inline auto operator-(T x, const fixed& y) noexcept -> fixed +{ + return fixed(x) -= y; +} + +// +// Multiplication +// + +template +constexpr inline auto operator*(const fixed& x, const fixed& y) noexcept -> fixed +{ + return fixed(x) *= y; +} + +template ::value>::type* = nullptr> +constexpr inline auto operator*(const fixed& x, T y) noexcept -> fixed +{ + return fixed(x) *= y; +} + +template ::value>::type* = nullptr> +constexpr inline auto operator*(T x, const fixed& y) noexcept -> fixed +{ + return fixed(y) *= x; +} + +// +// Division +// + +template +constexpr inline auto operator/(const fixed& x, const fixed& y) noexcept -> fixed +{ + return fixed(x) /= y; +} + +template ::value>::type* = nullptr> +constexpr inline auto operator/(const fixed& x, T y) noexcept -> fixed +{ + return fixed(x) /= y; +} + +template ::value>::type* = nullptr> +constexpr inline auto operator/(T x, const fixed& y) noexcept -> fixed +{ + return fixed(x) /= y; +} + +// +// Shift operators +// + +template ::value>::type* = nullptr> +constexpr inline auto operator>>(const fixed& x, T y) noexcept -> fixed +{ + return fixed(x) >>= y; +} + +template ::value>::type* = nullptr> +constexpr inline auto operator<<(const fixed& x, T y) noexcept -> fixed +{ + return fixed(x) <<= y; +} + +// +// Comparison operators +// + +template +constexpr inline auto operator==(const fixed& x, const fixed& y) noexcept -> bool +{ + return x.raw_value() == y.raw_value(); +} + +template +constexpr inline auto operator!=(const fixed& x, const fixed& y) noexcept -> bool +{ + return x.raw_value() != y.raw_value(); +} + +template +constexpr inline auto operator<(const fixed& x, const fixed& y) noexcept -> bool +{ + return x.raw_value() < y.raw_value(); +} + +template +constexpr inline auto operator>(const fixed& x, const fixed& y) noexcept -> bool +{ + return x.raw_value() > y.raw_value(); +} + +template +constexpr inline auto operator<=(const fixed& x, const fixed& y) noexcept -> bool +{ + return x.raw_value() <= y.raw_value(); +} + +template +constexpr inline auto operator>=(const fixed& x, const fixed& y) noexcept -> bool +{ + return x.raw_value() >= y.raw_value(); +} + +namespace detail +{ + +QT_WARNING_PUSH +QT_WARNING_DISABLE_CLANG("-Wunneeded-internal-declaration") + +// Number of base-10 digits required to fully represent a number of bits +static constexpr auto max_digits10(int bits) -> int +{ + // 8.24 fixed-point equivalent of (int)ceil(bits * std::log10(2)); + using T = long long; // NOLINT(google-runtime-int) + return static_cast((T{bits} * 5050445 + (T{1} << 24) - 1) >> 24); // NOLINT(hicpp-signed-bitwise) +} + +// Number of base-10 digits that can be fully represented by a number of bits +static constexpr auto digits10(int bits) -> int +{ + // 8.24 fixed-point equivalent of (int)(bits * std::log10(2)); + using T = long long; // NOLINT(google-runtime-int) + return static_cast((T{bits} * 5050445) >> 24); //NOLINT(hicpp-signed-bitwise) +} + +QT_WARNING_POP + +} // namespace detail +} // namespace fpm + +// Specializations for customization points +namespace std // NOLINT(cert-dcl58-cpp) +{ + +template +struct hash> +{ + using argument_type = fpm::fixed; + using result_type = std::size_t; + + auto operator()(const argument_type &arg) const noexcept(noexcept(std::declval>()(arg.raw_value()))) + -> result_type + { + return m_hash(arg.raw_value()); + } + +private: + std::hash m_hash; +}; + +template +struct numeric_limits> +{ + static constexpr bool is_specialized = true; + static constexpr bool is_signed = std::numeric_limits::is_signed; + static constexpr bool is_integer = false; + static constexpr bool is_exact = true; + static constexpr bool has_infinity = false; + static constexpr bool has_quiet_NaN = false; + static constexpr bool has_signaling_NaN = false; + static constexpr bool has_denorm = std::denorm_absent; + static constexpr bool has_denorm_loss = false; + static constexpr std::float_round_style round_style = std::round_to_nearest; + static constexpr bool is_iec_559 = false; + static constexpr bool is_bounded = true; + static constexpr bool is_modulo = std::numeric_limits::is_modulo; + static constexpr int digits = std::numeric_limits::digits; + + // Any number with `digits10` significant base-10 digits (that fits in + // the range of the type) is guaranteed to be convertible from text and + // back without change. Worst case, this is 0.000...001, so we can only + // guarantee this case. Nothing more. + static constexpr int digits10 = 1; + + // This is equal to max_digits10 for the integer and fractional part together. + static constexpr int max_digits10 = + fpm::detail::max_digits10(std::numeric_limits::digits - F) + fpm::detail::max_digits10(F); + + static constexpr int radix = 2; + static constexpr int min_exponent = 1 - F; + static constexpr int min_exponent10 = -fpm::detail::digits10(F); + static constexpr int max_exponent = std::numeric_limits::digits - F; + static constexpr int max_exponent10 = fpm::detail::digits10(std::numeric_limits::digits - F); + static constexpr bool traps = true; + static constexpr bool tinyness_before = false; + + static constexpr auto lowest() noexcept -> fpm::fixed + { + return fpm::fixed::from_raw_value(std::numeric_limits::lowest()); + } + + static constexpr auto min() noexcept -> fpm::fixed + { + return lowest(); + } + + static constexpr auto max() noexcept -> fpm::fixed + { + return fpm::fixed::from_raw_value(std::numeric_limits::max()); + } + + static constexpr auto epsilon() noexcept -> fpm::fixed + { + return fpm::fixed::from_raw_value(1); + } + + static constexpr auto round_error() noexcept -> fpm::fixed + { + return fpm::fixed(1) / 2; + } + + static constexpr auto denorm_min() noexcept -> fpm::fixed + { + return min(); + } +}; + +} // namespace std + +#endif diff --git a/src/libs/vmisc/fpm/math.hpp b/src/libs/vmisc/fpm/math.hpp new file mode 100644 index 000000000..d41110cc6 --- /dev/null +++ b/src/libs/vmisc/fpm/math.hpp @@ -0,0 +1,700 @@ +#ifndef FPM_MATH_HPP +#define FPM_MATH_HPP + +#include "fixed.hpp" +#include + +#ifdef _MSC_VER +#include +#endif + +#include + +#if QT_VERSION < QT_VERSION_CHECK(5, 5, 0) +#include "../diagnostic.h" +#endif // QT_VERSION < QT_VERSION_CHECK(5, 5, 0) + +namespace fpm +{ + +// +// Helper functions +// +namespace detail +{ + +QT_WARNING_PUSH +QT_WARNING_DISABLE_CLANG("-Wsign-conversion") + +// Returns the index of the most-signifcant set bit +inline auto find_highest_bit(unsigned long long value) noexcept -> long +{ + assert(value != 0); +#if defined(_MSC_VER) + unsigned long index; +#if defined(_WIN64) + _BitScanReverse64(&index, value); +#else + if (_BitScanReverse(&index, static_cast(value >> 32)) != 0) { + index += 32; + } else { + _BitScanReverse(&index, static_cast(value & 0xfffffffflu)); + } +#endif + return index; +#elif defined(__GNUC__) || defined(__clang__) + return sizeof(value) * 8 - 1 - __builtin_clzll(value); // NOLINT +#else +# error "your platform does not support find_highest_bit()" +#endif +} + +QT_WARNING_POP + +} // namespace detail + +// +// Classification methods +// + +template +constexpr inline auto fpclassify(fixed x) noexcept -> int +{ + return (x.raw_value() == 0) ? FP_ZERO : FP_NORMAL; +} + +template +constexpr inline auto isfinite(fixed /*unused*/) noexcept -> bool +{ + return true; +} + +template +constexpr inline auto isinf(fixed /*unused*/) noexcept -> bool +{ + return false; +} + +template +constexpr inline auto isnan(fixed /*unused*/) noexcept -> bool +{ + return false; +} + +template +constexpr inline auto isnormal(fixed x) noexcept -> bool +{ + return x.raw_value() != 0; +} + +template +constexpr inline auto signbit(fixed x) noexcept -> bool +{ + return x.raw_value() < 0; +} + +template +constexpr inline auto isgreater(fixed x, fixed y) noexcept -> bool +{ + return x > y; +} + +template +constexpr inline auto isgreaterequal(fixed x, fixed y) noexcept -> bool +{ + return x >= y; +} + +template +constexpr inline auto isless(fixed x, fixed y) noexcept -> bool +{ + return x < y; +} + +template +constexpr inline auto islessequal(fixed x, fixed y) noexcept -> bool +{ + return x <= y; +} + +template +constexpr inline auto islessgreater(fixed x, fixed y) noexcept -> bool +{ + return x != y; +} + +template +constexpr inline auto isunordered(fixed /*x*/, fixed /*y*/) noexcept -> bool +{ + return false; +} + +// +// Nearest integer operations +// +template +inline auto ceil(fixed x) noexcept -> fixed +{ + constexpr auto FRAC = B(1) << F; + auto value = x.raw_value(); + if (value > 0) + { + value += FRAC - 1; + } + return fixed::from_raw_value(value / FRAC * FRAC); +} + +template +inline auto floor(fixed x) noexcept -> fixed +{ + constexpr auto FRAC = B(1) << F; + auto value = x.raw_value(); + if (value < 0) + { + value -= FRAC - 1; + } + return fixed::from_raw_value(value / FRAC * FRAC); +} + +template +inline auto trunc(fixed x) noexcept -> fixed +{ + constexpr auto FRAC = B(1) << F; + return fixed::from_raw_value(x.raw_value() / FRAC * FRAC); +} + +template +inline auto round(fixed x) noexcept -> fixed +{ + constexpr auto FRAC = B(1) << F; + auto value = x.raw_value() / (FRAC / 2); + return fixed::from_raw_value(((value / 2) + (value % 2)) * FRAC); +} + +template +auto nearbyint(fixed x) noexcept -> fixed +{ + // Rounding mode is assumed to be FE_TONEAREST + constexpr auto FRAC = B(1) << F; + auto value = x.raw_value(); + const bool is_half = std::abs(value % FRAC) == FRAC / 2; + value /= FRAC / 2; + value = (value / 2) + (value % 2); + value -= (value % 2) * is_half; + return fixed::from_raw_value(value * FRAC); +} + +template +constexpr inline auto rint(fixed x) noexcept -> fixed +{ + // Rounding mode is assumed to be FE_TONEAREST + return nearbyint(x); +} + +// +// Mathematical functions +// +template +constexpr inline auto abs(fixed x) noexcept -> fixed +{ + return (x >= fixed{0}) ? x : -x; +} + +template +constexpr inline auto fmod(fixed x, fixed y) noexcept -> fixed +{ + return + assert(y.raw_value() != 0), + fixed::from_raw_value(x.raw_value() % y.raw_value()); +} + +template +constexpr inline auto remainder(fixed x, fixed y) noexcept -> fixed +{ + return + assert(y.raw_value() != 0), + x - nearbyint(x / y) * y; +} + +template +inline auto remquo(fixed x, fixed y, int* quo) noexcept -> fixed +{ + assert(y.raw_value() != 0); + assert(quo != nullptr); + *quo = x.raw_value() / y.raw_value(); + return fixed::from_raw_value(x.raw_value() % y.raw_value()); +} + +// +// Manipulation functions +// + +template +constexpr inline auto copysign(fixed x, fixed y) noexcept -> fixed +{ + return + x = abs(x), + (y >= fixed{0}) ? x : -x; +} + +template +constexpr inline auto nextafter(fixed from, fixed to) noexcept -> fixed +{ + return from == to ? to : + to > from ? fixed::from_raw_value(from.raw_value() + 1) + : fixed::from_raw_value(from.raw_value() - 1); +} + +template +constexpr inline auto nexttoward(fixed from, fixed to) noexcept -> fixed +{ + return nextafter(from, to); +} + +template +inline auto modf(fixed x, fixed* iptr) noexcept -> fixed +{ + const auto raw = x.raw_value(); + constexpr auto FRAC = B{1} << F; + *iptr = fixed::from_raw_value(raw / FRAC * FRAC); + return fixed::from_raw_value(raw % FRAC); +} + + +// +// Power functions +// + +template ::value>::type* = nullptr> +auto pow(fixed base, T exp) noexcept -> fixed +{ + using Fixed = fixed; + + if (base == Fixed(0)) { + assert(exp > 0); + return Fixed(0); + } + + Fixed result {1}; + if (exp < 0) + { + for (Fixed intermediate = base; exp != 0; exp /= 2, intermediate *= intermediate) + { + if ((exp % 2) != 0) + { + result /= intermediate; + } + } + } + else + { + for (Fixed intermediate = base; exp != 0; exp /= 2, intermediate *= intermediate) + { + if ((exp % 2) != 0) + { + result *= intermediate; + } + } + } + return result; +} + +template +auto pow(fixed base, fixed exp) noexcept -> fixed +{ + using Fixed = fixed; + + if (base == Fixed(0)) { + assert(exp > Fixed(0)); + return Fixed(0); + } + + if (exp < Fixed(0)) + { + return 1 / pow(base, -exp); + } + + constexpr auto FRAC = B(1) << F; + if (exp.raw_value() % FRAC == 0) + { + // Non-fractional exponents are easier to calculate + return pow(base, exp.raw_value() / FRAC); + } + + // For negative bases we do not support fractional exponents. + // Technically fractions with odd denominators could work, + // but that's too much work to figure out. + assert(base > Fixed(0)); + return exp2(log2(base) * exp); +} + +template +auto exp(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + if (x < Fixed(0)) { + return 1 / exp(-x); + } + constexpr auto FRAC = B(1) << F; + const B x_int = x.raw_value() / FRAC; + x -= x_int; + assert(x >= Fixed(0) && x < Fixed(1)); + + constexpr auto fA = Fixed::template from_fixed_point<63>( 128239257017632854LL); // 1.3903728105644451e-2 + constexpr auto fB = Fixed::template from_fixed_point<63>( 320978614890280666LL); // 3.4800571158543038e-2 + constexpr auto fC = Fixed::template from_fixed_point<63>(1571680799599592947LL); // 1.7040197373796334e-1 + constexpr auto fD = Fixed::template from_fixed_point<63>(4603349000587966862LL); // 4.9909609871464493e-1 + constexpr auto fE = Fixed::template from_fixed_point<62>(4612052447974689712LL); // 1.0000794567422495 + constexpr auto fF = Fixed::template from_fixed_point<63>(9223361618412247875LL); // 9.9999887043019773e-1 + return pow(Fixed::e(), x_int) * (((((fA * x + fB) * x + fC) * x + fD) * x + fE) * x + fF); +} + +template +auto exp2(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + if (x < Fixed(0)) { + return 1 / exp2(-x); + } + constexpr auto FRAC = B(1) << F; + const B x_int = x.raw_value() / FRAC; + x -= x_int; + assert(x >= Fixed(0) && x < Fixed(1)); + + constexpr auto fA = Fixed::template from_fixed_point<63>( 17491766697771214LL); // 1.8964611454333148e-3 + constexpr auto fB = Fixed::template from_fixed_point<63>( 82483038782406547LL); // 8.9428289841091295e-3 + constexpr auto fC = Fixed::template from_fixed_point<63>( 515275173969157690LL); // 5.5866246304520701e-2 + constexpr auto fD = Fixed::template from_fixed_point<63>(2214897896212987987LL); // 2.4013971109076949e-1 + constexpr auto fE = Fixed::template from_fixed_point<63>(6393224161192452326LL); // 6.9315475247516736e-1 + constexpr auto fF = Fixed::template from_fixed_point<63>(9223371050976163566LL); // 9.9999989311082668e-1 + return Fixed(1 << x_int) * (((((fA * x + fB) * x + fC) * x + fD) * x + fE) * x + fF); +} + +template +auto expm1(fixed x) noexcept -> fixed +{ + return exp(x) - 1; // cppcheck-suppress unpreciseMathCall +} + +template +auto log2(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + assert(x > Fixed(0)); + + // Normalize input to the [1:2] domain + B value = x.raw_value(); + const long highest = detail::find_highest_bit(value); + if (highest >= F) + { + value >>= (highest - F); + } + else + { + value <<= (F - highest); + } + x = Fixed::from_raw_value(value); + assert(x >= Fixed(1) && x < Fixed(2)); + + constexpr auto fA = Fixed::template from_fixed_point<63>( 413886001457275979LL); // 4.4873610194131727e-2 + constexpr auto fB = Fixed::template from_fixed_point<63>(-3842121857793256941LL); // -4.1656368651734915e-1 + constexpr auto fC = Fixed::template from_fixed_point<62>( 7522345947206307744LL); // 1.6311487636297217 + constexpr auto fD = Fixed::template from_fixed_point<61>(-8187571043052183818LL); // -3.5507929249026341 + constexpr auto fE = Fixed::template from_fixed_point<60>( 5870342889289496598LL); // 5.0917108110420042 + constexpr auto fF = Fixed::template from_fixed_point<61>(-6457199832668582866LL); // -2.8003640347009253 + return Fixed(highest - F) + (((((fA * x + fB) * x + fC) * x + fD) * x + fE) * x + fF); +} + +template +auto log(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + return log2(x) / log2(Fixed::e()); +} + +template +auto log10(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + return log2(x) / log2(Fixed(10)); +} + +template +auto log1p(fixed x) noexcept -> fixed +{ + return log(1 + x); // cppcheck-suppress unpreciseMathCall +} + +template +auto cbrt(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + + if (x == Fixed(0)) + { + return x; + } + if (x < Fixed(0)) + { + return -cbrt(-x); + } + assert(x >= Fixed(0)); + + // Finding the cube root of an integer, taken from Hacker's Delight, + // based on the square root algorithm. + + // We start at the greatest power of eight that's less than the argument. + int ofs = ((detail::find_highest_bit(x.raw_value()) + 2*F) / 3 * 3); + I num = I{x.raw_value()}; + I res = 0; + + const auto do_round = [&] + { + for (; ofs >= 0; ofs -= 3) + { + res += res; + const I val = (3*res*(res + 1) + 1) << ofs; + if (num >= val) + { + num -= val; + res++; + } + } + }; + + // We should shift by 2*F (since there are two multiplications), but that + // could overflow even the intermediate type, so we have to split the + // algorithm up in two rounds of F bits each. Each round will deplete + // 'num' digit by digit, so after a round we can shift it again. + num <<= F; + ofs -= F; + do_round(); + + num <<= F; + ofs += F; + do_round(); + + return Fixed::from_raw_value(static_cast(res)); +} + +template +auto sqrt(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + + assert(x >= Fixed(0)); + if (x == Fixed(0)) + { + return x; + } + + // Finding the square root of an integer in base-2, from: + // https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29 + + // Shift by F first because it's fixed-point. + I num = I{x.raw_value()} << F; + I res = 0; + + // "bit" starts at the greatest power of four that's less than the argument. + for (I bit = I{1} << ((detail::find_highest_bit(x.raw_value()) + F) / 2 * 2); bit != 0; bit >>= 2) + { + const I val = res + bit; + res >>= 1; + if (num >= val) + { + num -= val; + res += bit; + } + } + + // Round the last digit up if necessary + if (num > res) + { + res++; + } + + return Fixed::from_raw_value(static_cast(res)); +} + +template +auto hypot(fixed x, fixed y) noexcept -> fixed +{ + assert(x != 0 || y != 0); + return sqrt(x*x + y*y); +} + +// +// Trigonometry functions +// + +template +auto sin(fixed x) noexcept -> fixed +{ + // This sine uses a fifth-order curve-fitting approximation originally + // described by Jasper Vijn on coranac.com which has a worst-case + // relative error of 0.07% (over [-pi:pi]). + using Fixed = fixed; + + // Turn x from [0..2*PI] domain into [0..4] domain + x = fmod(x, Fixed::two_pi()); + x = x / Fixed::half_pi(); + + // Take x modulo one rotation, so [-4..+4]. + if (x < Fixed(0)) { + x += Fixed(4); + } + + int sign = +1; + if (x > Fixed(2)) { + // Reduce domain to [0..2]. + sign = -1; + x -= Fixed(2); + } + + if (x > Fixed(1)) { + // Reduce domain to [0..1]. + x = Fixed(2) - x; + } + + const Fixed x2 = x*x; + return sign * x * (Fixed::pi() - x2*(Fixed::two_pi() - 5 - x2*(Fixed::pi() - 3)))/2; +} + +template +inline auto cos(fixed x) noexcept -> fixed +{ + return sin(fixed::half_pi() + x); +} + +template +inline auto tan(fixed x) noexcept -> fixed +{ + auto cx = cos(x); + + // Tangent goes to infinity at 90 and -90 degrees. + // We can't represent that with fixed-point maths. + assert(abs(cx).raw_value() > 1); + + return sin(x) / cx; +} + +namespace detail { + +// Calculates atan(x) assuming that x is in the range [0,1] +template +auto atan_sanitized(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + assert(x >= Fixed(0) && x <= Fixed(1)); + + constexpr auto fA = Fixed::template from_fixed_point<63>( 716203666280654660LL); // 0.0776509570923569 + constexpr auto fB = Fixed::template from_fixed_point<63>(-2651115102768076601LL); // -0.287434475393028 + constexpr auto fC = Fixed::template from_fixed_point<63>( 9178930894564541004LL); // 0.995181681698119 (PI/4 - A - B) + + const auto xx = x * x; + return ((fA*xx + fB)*xx + fC)*x; +} + +// Calculate atan(y / x), assuming x != 0. +// +// If x is very, very small, y/x can easily overflow the fixed-point range. +// If q = y/x and q > 1, atan(q) would calculate atan(1/q) as intermediate step +// anyway. We can shortcut that here and avoid the loss of information, thus +// improving the accuracy of atan(y/x) for very small x. +template +auto atan_div(fixed y, fixed x) noexcept -> fixed +{ + using Fixed = fixed; + assert(x != Fixed(0)); + + // Make sure y and x are positive. + // If y / x is negative (when y or x, but not both, are negative), negate the result to + // keep the correct outcome. + if (y < Fixed(0)) { + if (x < Fixed(0)) { + return atan_div(-y, -x); + } + return -atan_div(-y, x); + } + if (x < Fixed(0)) { + return -atan_div(y, -x); + } + assert(y >= Fixed(0)); + assert(x > Fixed(0)); + + if (y > x) { + return Fixed::half_pi() - detail::atan_sanitized(x / y); + } + return detail::atan_sanitized(y / x); +} + +} // namespace detail + +template +auto atan(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + if (x < Fixed(0)) + { + return -atan(-x); + } + + if (x > Fixed(1)) + { + return Fixed::half_pi() - detail::atan_sanitized(Fixed(1) / x); + } + + return detail::atan_sanitized(x); +} + +template +auto asin(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + assert(x >= Fixed(-1) && x <= Fixed(+1)); + + const auto yy = Fixed(1) - x * x; + if (yy == Fixed(0)) + { + return copysign(Fixed::half_pi(), x); + } + return detail::atan_div(x, sqrt(yy)); +} + +template +auto acos(fixed x) noexcept -> fixed +{ + using Fixed = fixed; + assert(x >= Fixed(-1) && x <= Fixed(+1)); + + if (x == Fixed(-1)) + { + return Fixed::pi(); + } + const auto yy = Fixed(1) - x * x; + return Fixed(2)*detail::atan_div(sqrt(yy), Fixed(1) + x); +} + +template +auto atan2(fixed y, fixed x) noexcept -> fixed +{ + using Fixed = fixed; + if (x == Fixed(0)) + { + assert(y != Fixed(0)); + return (y > Fixed(0)) ? Fixed::half_pi() : -Fixed::half_pi(); + } + + auto ret = detail::atan_div(y, x); + + if (x < Fixed(0)) + { + return (y >= Fixed(0)) ? ret + Fixed::pi() : ret - Fixed::pi(); + } + return ret; +} + +} // namespace fpm + +#endif diff --git a/src/libs/vmisc/vmisc.pri b/src/libs/vmisc/vmisc.pri index c8541cd0e..2c1e44094 100644 --- a/src/libs/vmisc/vmisc.pri +++ b/src/libs/vmisc/vmisc.pri @@ -55,7 +55,9 @@ HEADERS += \ $$PWD/vmodifierkey.h \ $$PWD/typedef.h \ $$PWD/backport/qscopeguard.h \ - $$PWD/dialogs/dialogselectlanguage.h + $$PWD/dialogs/dialogselectlanguage.h \ + $$PWD/fpm/fixed.hpp \ + $$PWD/fpm/math.hpp contains(DEFINES, APPIMAGE) { SOURCES += \ diff --git a/src/libs/vtools/dialogs/tools/dialogsinglepoint.cpp b/src/libs/vtools/dialogs/tools/dialogsinglepoint.cpp index b9ec1ced9..649ce8464 100644 --- a/src/libs/vtools/dialogs/tools/dialogsinglepoint.cpp +++ b/src/libs/vtools/dialogs/tools/dialogsinglepoint.cpp @@ -53,8 +53,8 @@ DialogSinglePoint::DialogSinglePoint(const VContainer *data, quint32 toolId, QWi ui->lineEditName->setClearButtonEnabled(true); - ui->doubleSpinBoxX->setRange(0, VAbstractValApplication::VApp()->fromPixel(SceneSize)); - ui->doubleSpinBoxY->setRange(0, VAbstractValApplication::VApp()->fromPixel(SceneSize)); + ui->doubleSpinBoxX->setRange(0, VAbstractValApplication::VApp()->fromPixel(maxSceneSize)); + ui->doubleSpinBoxY->setRange(0, VAbstractValApplication::VApp()->fromPixel(maxSceneSize)); InitOkCancel(ui); connect(ui->lineEditName, &QLineEdit::textChanged, this, [this]() diff --git a/src/libs/vwidgets/vmaingraphicsview.cpp b/src/libs/vwidgets/vmaingraphicsview.cpp index 5ae5a1415..f33aff0bb 100644 --- a/src/libs/vwidgets/vmaingraphicsview.cpp +++ b/src/libs/vwidgets/vmaingraphicsview.cpp @@ -51,21 +51,18 @@ #include #include #include +#include #include "../vmisc/def.h" -#include "../vmisc/vmath.h" #include "vmaingraphicsscene.h" #include "vsimplecurve.h" #include "vcontrolpointspline.h" #include "../vmisc/vabstractapplication.h" #include "../vmisc/vcommonsettings.h" #include "../vmisc/literals.h" -#include "vabstractmainwindow.h" #include "global.h" #include "../ifc/xml/utils.h" -const qreal maxSceneSize = ((20.0 * 1000.0) / 25.4) * PrintDPI; // 20 meters in pixels - namespace { qreal ScrollingSteps(QWheelEvent* wheel_event)