1 /*
2     SPDX-FileCopyrightText: 2010 Johannes Loehnert <loehnert.kde@gmx.de>
3 
4     SPDX-License-Identifier: GPL-2.0-or-later
5 */
6 
7 #include "utilities.h"
8 
9 #include <QDebug>
10 #include <QRandomGenerator>
11 
getBestFit(int & xCount,int & yCount,qreal target_aspect,int approx_count)12 void getBestFit(int &xCount, int &yCount, qreal target_aspect, int approx_count) {
13     qreal nx_exact = sqrt(approx_count * target_aspect);
14     // avoid taking the sqrt again
15     qreal ny_exact = approx_count / nx_exact;
16 
17     // catch very odd cases
18     if (nx_exact < 1) nx_exact = 1.01;
19     if (ny_exact < 1) ny_exact = 1.01;
20 
21     qreal aspect1 = floor(nx_exact) / ceil(ny_exact);
22     qreal aspect2 = ceil(nx_exact) / floor(ny_exact);
23 
24     aspect1 = target_aspect - aspect1;
25     aspect2 = aspect2 - target_aspect;
26 
27     if (aspect1 < aspect2) ny_exact += 1.0; else nx_exact += 1.0;
28 
29     xCount = (int)(floor(nx_exact) + 0.1);
30     yCount = (int)(floor(ny_exact) + 0.1);
31 }
32 
getBestFitExtended(int & xCount,int & yCount,qreal target_aspect,int approx_count,qreal tiles_per_cell,qreal additional_tiles_per_row,qreal additional_tiles_per_column,qreal additional_tiles)33 void getBestFitExtended(int &xCount, int &yCount, qreal target_aspect, int approx_count,
34             qreal tiles_per_cell, qreal additional_tiles_per_row, qreal additional_tiles_per_column, qreal additional_tiles) {
35     // solves the equations
36     //  N = TPC * x * y  +  ATPC * x + ATPR * y + AT
37     //  target_aspect = x / y
38     //
39     // for x, y; and rounds them to the nearest integer values giving least distance to target_aspect.
40     qreal p_half = (target_aspect * additional_tiles_per_column + additional_tiles_per_row) / (2 * target_aspect * tiles_per_cell);
41     qreal q = (approx_count - additional_tiles) / (target_aspect * tiles_per_cell);
42 
43     qreal p_half_sq = p_half*p_half;
44     if (p_half_sq + q < 0) {
45         xCount = 1;
46         yCount = 1;
47         return;
48     }
49 
50     qreal ny_exact = -p_half + sqrt(p_half_sq + q);
51     qreal nx_exact = target_aspect * ny_exact;
52 
53     qDebug() << "nx_exact: " << nx_exact << " ny_exact: " << ny_exact <<
54                 "giving count: " << (tiles_per_cell*nx_exact*ny_exact + additional_tiles_per_column * nx_exact + additional_tiles_per_row * ny_exact + additional_tiles);
55 
56     // catch very odd cases
57     if (nx_exact < 1) nx_exact = 1.01;
58     if (ny_exact < 1) ny_exact = 1.01;
59 
60     qreal aspect1 = floor(nx_exact) / ceil(ny_exact);
61     qreal aspect2 = ceil(nx_exact) / floor(ny_exact);
62     qreal aspect3 = ceil(nx_exact) / ceil(ny_exact);
63 
64     aspect1 = target_aspect - aspect1;
65     aspect2 = aspect2 - target_aspect;
66     aspect3 = abs(aspect3 - target_aspect);
67 
68     if (aspect1 < aspect2) {
69         ny_exact += 1.0;
70         if (aspect3 < aspect1) nx_exact += 1.0;
71     }
72     else {
73         nx_exact += 1.0;
74         if (aspect3 < aspect2) ny_exact += 1.0;
75     }
76 
77     xCount = (int)(floor(nx_exact) + 0.1);
78     yCount = (int)(floor(ny_exact) + 0.1);
79 }
80 
81 // skews x with "strength" a.
82 // x is expected to lie within [0, 1].
83 // a = +/-1 is already a very strong skew.
84 // negative a: skew towards x=0, positive: skew towards x=1.
skew_randnum(qreal x,qreal a)85 qreal skew_randnum(qreal x, qreal a) {
86     if (a==0) return x;
87 
88     qreal asq = exp(-2 * abs(a));
89     if (a>0) x = 1-x;
90 
91     qreal mp2 = (x-1) * (2/asq - 1);
92     qreal q = (x-1)*(x-1) - 1;
93 
94     // We apply a function on x, which is a hyperbola through (0,0) and (1,1)
95     // with (1, 0) as focal point. You really don't want to know the gory details.
96     x = mp2 + sqrt(mp2*mp2 - q);
97 
98     if (a>0) x = 1-x;
99     return x;
100 }
101 
102 
nonuniform_rand(qreal min,qreal max,qreal sigma,qreal skew)103 qreal nonuniform_rand(qreal min, qreal max, qreal sigma, qreal skew) {
104 
105     // 0.4247: sigma at which distribution function is 1/2 of max at interval boundaries
106 
107     qreal randNum;
108 
109     auto *generator = QRandomGenerator::global();
110     if (sigma > 0.4247) {
111         // "wide" distribution, use rejection sampling
112         qreal x, y;
113         qreal ssq = 2 * sigma * sigma;
114         do {
115             x = 0.000001 * qreal(generator->bounded(1000000));
116             y = 0.000001 * qreal(generator->bounded(1000000));
117         } while (y > exp(-(x-0.5)*(x-0.5)/ssq));
118 
119         randNum = x;
120     }
121     else {
122         // "narrow" distribution, use Marsaglia method until a random number within 0, 1 is found.
123         qreal u1, u2, p, q, x1, x2;
124 
125         randNum = -1;
126         do {
127             do {
128                 u1 = 0.000002 * qreal(generator->bounded(1000000)) - 1;
129                 u2 = 0.000002 * qreal(generator->bounded(1000000)) - 1;
130                 q = u1*u1 + u2*u2;
131             } while (q>1);
132             p = sqrt(-2 * log(q) / q) * sigma;
133             x1 = u1 * p + 0.5;
134             x2 = u2 * p + 0.5;
135 
136             if (x1>= 0 && x1 <= 1) {
137                 randNum = x1;
138             }
139             else {
140                 if (x2 >= 0 && x2 <= 1) randNum = x2;
141             }
142         } while (randNum < 0);
143     }
144 
145     return min + (max - min) * skew_randnum(randNum, skew);
146 }
147 
148