MRTKLIB 0.7.14
Modernized GNSS Positioning Library
Loading...
Searching...
No Matches
mrtk_lambda.c
Go to the documentation of this file.
1/*------------------------------------------------------------------------------
2 * mrtk_lambda.c : integer ambiguity resolution (LAMBDA/MLAMBDA)
3 *
4 * Copyright (C) 2026 H.SHIONO (MRTKLIB Project)
5 * Copyright (C) 2023-2025 Cabinet Office, Japan
6 * Copyright (C) 2024-2025 Lighthouse Technology & Consulting Co. Ltd.
7 * Copyright (C) 2023-2025 Japan Aerospace Exploration Agency
8 * Copyright (C) 2023-2025 TOSHIBA ELECTRONIC TECHNOLOGIES CORPORATION
9 * Copyright (C) 2015- Mitsubishi Electric Corp.
10 * Copyright (C) 2014 Geospatial Information Authority of Japan
11 * Copyright (C) 2014 T.SUZUKI
12 * Copyright (C) 2007-2023 T.TAKASU
13 *
14 * SPDX-License-Identifier: BSD-2-Clause
15 *----------------------------------------------------------------------------*/
16#include "mrtklib/mrtk_lambda.h"
17
18#include <math.h>
19#include <stdio.h>
20#include <stdlib.h>
21#include <string.h>
22
23#include "mrtklib/mrtk_mat.h"
24
25/* constants/macros ----------------------------------------------------------*/
26
27#define LOOPMAX 10000 /* maximum count of search loop */
28
29#define SGN(x) ((x) <= 0.0 ? -1.0 : 1.0)
30#define ROUND(x) (floor((x) + 0.5))
31#define SWAP(x, y) \
32 do { \
33 double tmp_; \
34 tmp_ = x; \
35 x = y; \
36 y = tmp_; \
37 } while (0)
38
39/* LD factorization (Q=L'*diag(D)*L) -----------------------------------------*/
40static int LD(int n, const double* Q, double* L, double* D) {
41 int i, j, k, info = 0;
42 double a, *A = mat(n, n);
43
44 memcpy(A, Q, sizeof(double) * n * n);
45 for (i = n - 1; i >= 0; i--) {
46 if ((D[i] = A[i + i * n]) <= 0.0) {
47 info = -1;
48 break;
49 }
50 a = sqrt(D[i]);
51 for (j = 0; j <= i; j++) {
52 L[i + j * n] = A[i + j * n] / a;
53 }
54 for (j = 0; j <= i - 1; j++) {
55 for (k = 0; k <= j; k++) {
56 A[j + k * n] -= L[i + k * n] * L[i + j * n];
57 }
58 }
59 for (j = 0; j <= i; j++) {
60 L[i + j * n] /= L[i + i * n];
61 }
62 }
63 free(A);
64 if (info) {
65 fprintf(stderr, "%s : LD factorization error\n", __FILE__);
66 }
67 return info;
68}
69/* integer gauss transformation ----------------------------------------------*/
70static void gauss(int n, double* L, double* Z, int i, int j) {
71 int k, mu;
72
73 if ((mu = (int)ROUND(L[i + j * n])) != 0) {
74 for (k = i; k < n; k++) {
75 L[k + n * j] -= (double)mu * L[k + i * n];
76 }
77 for (k = 0; k < n; k++) {
78 Z[k + n * j] -= (double)mu * Z[k + i * n];
79 }
80 }
81}
82/* permutations --------------------------------------------------------------*/
83static void perm(int n, double* L, double* D, int j, double del, double* Z) {
84 int k;
85 double eta, lam, a0, a1;
86
87 eta = D[j] / del;
88 lam = D[j + 1] * L[j + 1 + j * n] / del;
89 D[j] = eta * D[j + 1];
90 D[j + 1] = del;
91 for (k = 0; k <= j - 1; k++) {
92 a0 = L[j + k * n];
93 a1 = L[j + 1 + k * n];
94 L[j + k * n] = -L[j + 1 + j * n] * a0 + a1;
95 L[j + 1 + k * n] = eta * a0 + lam * a1;
96 }
97 L[j + 1 + j * n] = lam;
98 for (k = j + 2; k < n; k++) {
99 SWAP(L[k + j * n], L[k + (j + 1) * n]);
100 }
101 for (k = 0; k < n; k++) {
102 SWAP(Z[k + j * n], Z[k + (j + 1) * n]);
103 }
104}
105/* lambda reduction (z=Z'*a, Qz=Z'*Q*Z=L'*diag(D)*L) (ref.[1]) ---------------*/
106static void reduction(int n, double* L, double* D, double* Z) {
107 int i, j, k;
108 double del;
109
110 j = n - 2;
111 k = n - 2;
112 while (j >= 0) {
113 if (j <= k) {
114 for (i = j + 1; i < n; i++) {
115 gauss(n, L, Z, i, j);
116 }
117 }
118 del = D[j] + L[j + 1 + j * n] * L[j + 1 + j * n] * D[j + 1];
119 if (del + 1E-6 < D[j + 1]) { /* compared considering numerical error */
120 perm(n, L, D, j, del, Z);
121 k = j;
122 j = n - 2;
123 } else {
124 j--;
125 }
126 }
127}
128/* modified lambda (mlambda) search (ref. [2]) -------------------------------*/
129static int search(int n, int m, const double* L, const double* D, const double* zs, double* zn, double* s) {
130 int i, j, k, c, nn = 0, imax = 0;
131 double newdist, maxdist = 1E99, y;
132 double *S = zeros(n, n), *dist = mat(n, 1), *zb = mat(n, 1), *z = mat(n, 1), *step = mat(n, 1);
133
134 k = n - 1;
135 dist[k] = 0.0;
136 zb[k] = zs[k];
137 z[k] = ROUND(zb[k]);
138 y = zb[k] - z[k];
139 step[k] = SGN(y);
140 for (c = 0; c < LOOPMAX; c++) {
141 newdist = dist[k] + y * y / D[k];
142 if (newdist < maxdist) {
143 if (k != 0) {
144 dist[--k] = newdist;
145 for (i = 0; i <= k; i++) {
146 S[k + i * n] = S[k + 1 + i * n] + (z[k + 1] - zb[k + 1]) * L[k + 1 + i * n];
147 }
148 zb[k] = zs[k] + S[k + k * n];
149 z[k] = ROUND(zb[k]);
150 y = zb[k] - z[k];
151 step[k] = SGN(y);
152 } else {
153 if (nn < m) {
154 if (nn == 0 || newdist > s[imax]) {
155 imax = nn;
156 }
157 for (i = 0; i < n; i++) {
158 zn[i + nn * n] = z[i];
159 }
160 s[nn++] = newdist;
161 } else {
162 if (newdist < s[imax]) {
163 for (i = 0; i < n; i++) {
164 zn[i + imax * n] = z[i];
165 }
166 s[imax] = newdist;
167 for (i = imax = 0; i < m; i++) {
168 if (s[imax] < s[i]) {
169 imax = i;
170 }
171 }
172 }
173 maxdist = s[imax];
174 }
175 z[0] += step[0];
176 y = zb[0] - z[0];
177 step[0] = -step[0] - SGN(step[0]);
178 }
179 } else {
180 if (k == n - 1) {
181 break;
182 } else {
183 k++;
184 z[k] += step[k];
185 y = zb[k] - z[k];
186 step[k] = -step[k] - SGN(step[k]);
187 }
188 }
189 }
190 for (i = 0; i < m - 1; i++) { /* sort by s */
191 for (j = i + 1; j < m; j++) {
192 if (s[i] < s[j]) {
193 continue;
194 }
195 SWAP(s[i], s[j]);
196 for (k = 0; k < n; k++) {
197 SWAP(zn[k + i * n], zn[k + j * n]);
198 }
199 }
200 }
201 free(S);
202 free(dist);
203 free(zb);
204 free(z);
205 free(step);
206
207 if (c >= LOOPMAX) {
208 fprintf(stderr, "%s : search loop count overflow\n", __FILE__);
209 return -1;
210 }
211 return 0;
212}
213/* lambda/mlambda integer least-square estimation ------------------------------
214 * integer least-square estimation. reduction is performed by lambda (ref.[1]),
215 * and search by mlambda (ref.[2]).
216 * args : int n I number of float parameters
217 * int m I number of fixed solutions
218 * double *a I float parameters (n x 1)
219 * double *Q I covariance matrix of float parameters (n x n)
220 * double *F O fixed solutions (n x m)
221 * double *s O sum of squared residulas of fixed solutions (1 x m)
222 * return : status (0:ok,other:error)
223 * notes : matrix stored by column-major order (fortran convension)
224 *-----------------------------------------------------------------------------*/
225extern int lambda(int n, int m, const double* a, const double* Q, double* F, double* s) {
226 int info;
227 double *L, *D, *Z, *z, *E;
228
229 if (n <= 0 || m <= 0) {
230 return -1;
231 }
232 L = zeros(n, n);
233 D = mat(n, 1);
234 Z = eye(n);
235 z = mat(n, 1);
236 E = mat(n, m);
237
238 /* LD factorization */
239 if (!(info = LD(n, Q, L, D))) {
240 /* lambda reduction */
241 reduction(n, L, D, Z);
242 matmul("TN", n, 1, n, 1.0, Z, a, 0.0, z); /* z=Z'*a */
243
244 /* mlambda search */
245 if (!(info = search(n, m, L, D, z, E, s))) {
246 info = solve("T", Z, E, n, m, F); /* F=Z'\E */
247 }
248 }
249 free(L);
250 free(D);
251 free(Z);
252 free(z);
253 free(E);
254 return info;
255}
256/* lambda reduction ------------------------------------------------------------
257 * reduction by lambda (ref [1]) for integer least square
258 * args : int n I number of float parameters
259 * double *Q I covariance matrix of float parameters (n x n)
260 * double *Z O lambda reduction matrix (n x n)
261 * return : status (0:ok,other:error)
262 *-----------------------------------------------------------------------------*/
263extern int lambda_reduction(int n, const double* Q, double* Z) {
264 double *L, *D;
265 int i, j, info;
266
267 if (n <= 0) {
268 return -1;
269 }
270
271 L = zeros(n, n);
272 D = mat(n, 1);
273
274 for (i = 0; i < n; i++) {
275 for (j = 0; j < n; j++) {
276 Z[i + j * n] = i == j ? 1.0 : 0.0;
277 }
278 }
279 /* LD factorization */
280 if ((info = LD(n, Q, L, D))) {
281 free(L);
282 free(D);
283 return info;
284 }
285 /* lambda reduction */
286 reduction(n, L, D, Z);
287
288 free(L);
289 free(D);
290 return 0;
291}
292/* mlambda search --------------------------------------------------------------
293 * search by mlambda (ref [2]) for integer least square
294 * args : int n I number of float parameters
295 * int m I number of fixed solutions
296 * double *a I float parameters (n x 1)
297 * double *Q I covariance matrix of float parameters (n x n)
298 * double *F O fixed solutions (n x m)
299 * double *s O sum of squared residulas of fixed solutions (1 x m)
300 * return : status (0:ok,other:error)
301 *-----------------------------------------------------------------------------*/
302extern int lambda_search(int n, int m, const double* a, const double* Q, double* F, double* s) {
303 double *L, *D;
304 int info;
305
306 if (n <= 0 || m <= 0) {
307 return -1;
308 }
309
310 L = zeros(n, n);
311 D = mat(n, 1);
312
313 /* LD factorization */
314 if ((info = LD(n, Q, L, D))) {
315 free(L);
316 free(D);
317 return info;
318 }
319 /* mlambda search */
320 info = search(n, m, L, D, a, F, s);
321
322 free(L);
323 free(D);
324 return info;
325}
326/* ratio test threshold tables -----------------------------------------------*/
327static const float inv_ratio_tbl_1[] = {
328 0.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
329 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
330 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
331 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
332 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.0010, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
333 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
334 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
335 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
336 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
337 1.0000, 1.0000, 0.0012, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9600, 0.9600, 0.9600,
338 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9700,
339 0.9700, 0.9700, 0.9800, 0.9800, 0.9800, 0.9800, 0.9800, 0.9800, 0.9800, 0.9800, 0.9800, 0.9800, 0.9800, 0.9800,
340 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900,
341 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.9900, 0.0015, 0.8700, 0.8700, 0.8800,
342 0.8800, 0.8900, 0.8900, 0.8900, 0.8900, 0.9100, 0.9200, 0.9200, 0.9200, 0.9200, 0.9300, 0.9300, 0.9300, 0.9300,
343 0.9300, 0.9300, 0.9300, 0.9300, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9500, 0.9500, 0.9500, 0.9500,
344 0.9500, 0.9500, 0.9500, 0.9500, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9600,
345 0.9600, 0.9600, 0.9600, 0.9600, 0.9600, 0.9700, 0.9700, 0.9700, 0.9700, 0.9700, 0.9700, 0.9700, 0.9700, 0.9700,
346 0.9700, 0.9700, 0.9700, 0.9700, 0.0020, 0.7800, 0.7800, 0.8000, 0.8000, 0.8100, 0.8200, 0.8300, 0.8400, 0.8400,
347 0.8600, 0.8600, 0.8600, 0.8700, 0.8800, 0.8800, 0.8800, 0.8900, 0.8900, 0.9000, 0.9000, 0.9000, 0.9000, 0.9000,
348 0.9000, 0.9100, 0.9100, 0.9100, 0.9200, 0.9200, 0.9200, 0.9200, 0.9300, 0.9300, 0.9300, 0.9400, 0.9400, 0.9400,
349 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9400, 0.9500,
350 0.9500, 0.9500, 0.9500, 0.9500, 0.9500, 0.9500, 0.9500, 0.9500, 0.9500, 0.9500, 0.9500, 0.9600, 0.0050, 0.5400,
351 0.5400, 0.5700, 0.5700, 0.5900, 0.6400, 0.6400, 0.6800, 0.6800, 0.6900, 0.7100, 0.7200, 0.7300, 0.7400, 0.7500,
352 0.7500, 0.7700, 0.7800, 0.8000, 0.8000, 0.8000, 0.8000, 0.8000, 0.8100, 0.8200, 0.8200, 0.8200, 0.8400, 0.8400,
353 0.8400, 0.8400, 0.8500, 0.8500, 0.8600, 0.8600, 0.8600, 0.8600, 0.8700, 0.8800, 0.8800, 0.8800, 0.8800, 0.8800,
354 0.8800, 0.8800, 0.8800, 0.8800, 0.8900, 0.8900, 0.8900, 0.9100, 0.9100, 0.9100, 0.9100, 0.9100, 0.9100, 0.9100,
355 0.9100, 0.9100, 0.9100, 0.9100, 0.9100, 0.9100, 0.0100, 0.3600, 0.3800, 0.4100, 0.4300, 0.4500, 0.5100, 0.5200,
356 0.5700, 0.5700, 0.5900, 0.6100, 0.6200, 0.6400, 0.6500, 0.6700, 0.6700, 0.6800, 0.6900, 0.7200, 0.7200, 0.7200,
357 0.7200, 0.7200, 0.7400, 0.7600, 0.7600, 0.7700, 0.7900, 0.7900, 0.7900, 0.7900, 0.8100, 0.8100, 0.8200, 0.8200,
358 0.8200, 0.8200, 0.8300, 0.8300, 0.8400, 0.8400, 0.8500, 0.8500, 0.8500, 0.8500, 0.8500, 0.8500, 0.8500, 0.8500,
359 0.8500, 0.8700, 0.8800, 0.8800, 0.8800, 0.8800, 0.8800, 0.8800, 0.8800, 0.8800, 0.8800, 0.8800, 0.8800, 0.8800,
360 0.0150, 0.2700, 0.2900, 0.3200, 0.3600, 0.3800, 0.4300, 0.4600, 0.5100, 0.5200, 0.5300, 0.5500, 0.5600, 0.5800,
361 0.6000, 0.6200, 0.6200, 0.6400, 0.6500, 0.6700, 0.6700, 0.6900, 0.6900, 0.6900, 0.7100, 0.7400, 0.7400, 0.7400,
362 0.7600, 0.7600, 0.7600, 0.7700, 0.7800, 0.7800, 0.7900, 0.7900, 0.8000, 0.8000, 0.8100, 0.8100, 0.8200, 0.8200,
363 0.8300, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400, 0.8400, 0.8400, 0.8400, 0.8500, 0.8500, 0.8600, 0.8600, 0.8600,
364 0.8600, 0.8600, 0.8600, 0.8600, 0.8700, 0.8700, 0.8700, 0.8800, 0.0200, 0.2100, 0.2400, 0.2700, 0.3000, 0.3300,
365 0.3800, 0.4200, 0.4600, 0.4800, 0.4900, 0.5200, 0.5200, 0.5500, 0.5700, 0.5900, 0.5900, 0.6200, 0.6200, 0.6500,
366 0.6500, 0.6700, 0.6700, 0.6700, 0.7000, 0.7200, 0.7200, 0.7200, 0.7400, 0.7500, 0.7500, 0.7500, 0.7700, 0.7700,
367 0.7800, 0.7800, 0.7800, 0.7900, 0.7900, 0.8000, 0.8100, 0.8100, 0.8100, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300,
368 0.8300, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400, 0.8500, 0.8500, 0.8500, 0.8500, 0.8500, 0.8500, 0.8600, 0.8600,
369 0.8600, 0.8700, 0.0250, 0.1900, 0.2200, 0.2400, 0.2700, 0.2900, 0.3500, 0.3900, 0.4200, 0.4500, 0.4600, 0.4900,
370 0.4900, 0.5300, 0.5400, 0.5600, 0.5600, 0.6000, 0.6000, 0.6300, 0.6400, 0.6500, 0.6600, 0.6600, 0.6800, 0.7100,
371 0.7100, 0.7100, 0.7200, 0.7200, 0.7300, 0.7400, 0.7500, 0.7500, 0.7700, 0.7700, 0.7700, 0.7800, 0.7800, 0.7900,
372 0.7900, 0.8000, 0.8000, 0.8000, 0.8100, 0.8100, 0.8100, 0.8200, 0.8200, 0.8300, 0.8300, 0.8300, 0.8300, 0.8300,
373 0.8400, 0.8400, 0.8500, 0.8500, 0.8500, 0.8500, 0.8600, 0.8600, 0.8600, 0.8700, 0.0300, 0.1700, 0.2000, 0.2200,
374 0.2400, 0.2600, 0.3200, 0.3500, 0.4000, 0.4200, 0.4400, 0.4700, 0.4700, 0.5100, 0.5300, 0.5400, 0.5500, 0.5800,
375 0.5900, 0.6200, 0.6200, 0.6400, 0.6500, 0.6500, 0.6700, 0.7000, 0.7000, 0.7000, 0.7100, 0.7100, 0.7200, 0.7300,
376 0.7400, 0.7500, 0.7600, 0.7600, 0.7600, 0.7700, 0.7700, 0.7800, 0.7800, 0.7900, 0.7900, 0.8000, 0.8000, 0.8100,
377 0.8100, 0.8100, 0.8200, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400, 0.8500, 0.8500, 0.8500,
378 0.8500, 0.8600, 0.8600, 0.8600, 0.0350, 0.1600, 0.1800, 0.2000, 0.2200, 0.2400, 0.3000, 0.3300, 0.3800, 0.4000,
379 0.4200, 0.4500, 0.4500, 0.4900, 0.5100, 0.5300, 0.5400, 0.5700, 0.5800, 0.6100, 0.6100, 0.6300, 0.6400, 0.6500,
380 0.6600, 0.6900, 0.6900, 0.6900, 0.7100, 0.7100, 0.7200, 0.7200, 0.7400, 0.7400, 0.7500, 0.7500, 0.7600, 0.7600,
381 0.7700, 0.7700, 0.7800, 0.7800, 0.7900, 0.7900, 0.7900, 0.8000, 0.8100, 0.8100, 0.8100, 0.8100, 0.8200, 0.8200,
382 0.8200, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400, 0.8500, 0.8500, 0.8500, 0.8500, 0.8600, 0.8600, 0.0400, 0.1500,
383 0.1700, 0.1800, 0.2100, 0.2200, 0.2700, 0.3100, 0.3600, 0.3900, 0.4000, 0.4400, 0.4400, 0.4800, 0.5000, 0.5200,
384 0.5300, 0.5600, 0.5600, 0.6000, 0.6100, 0.6200, 0.6400, 0.6500, 0.6600, 0.6800, 0.6900, 0.6900, 0.7000, 0.7100,
385 0.7100, 0.7100, 0.7300, 0.7400, 0.7400, 0.7500, 0.7500, 0.7600, 0.7600, 0.7700, 0.7700, 0.7800, 0.7800, 0.7900,
386 0.7900, 0.7900, 0.8000, 0.8100, 0.8100, 0.8100, 0.8100, 0.8200, 0.8200, 0.8200, 0.8300, 0.8400, 0.8400, 0.8400,
387 0.8500, 0.8500, 0.8500, 0.8500, 0.8500, 0.8600, 0.0450, 0.1400, 0.1600, 0.1700, 0.1900, 0.2100, 0.2600, 0.2900,
388 0.3400, 0.3800, 0.3900, 0.4300, 0.4400, 0.4800, 0.4900, 0.5100, 0.5200, 0.5600, 0.5600, 0.5900, 0.6000, 0.6100,
389 0.6300, 0.6500, 0.6500, 0.6700, 0.6800, 0.6800, 0.7000, 0.7100, 0.7100, 0.7100, 0.7300, 0.7300, 0.7400, 0.7400,
390 0.7500, 0.7500, 0.7600, 0.7700, 0.7700, 0.7800, 0.7800, 0.7800, 0.7900, 0.7900, 0.7900, 0.8000, 0.8000, 0.8100,
391 0.8100, 0.8100, 0.8200, 0.8200, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400, 0.8400, 0.8400, 0.8500, 0.8500, 0.8600,
392 0.0500, 0.1300, 0.1500, 0.1600, 0.1800, 0.1900, 0.2400, 0.2800, 0.3200, 0.3700, 0.3900, 0.4200, 0.4300, 0.4700,
393 0.4800, 0.5100, 0.5200, 0.5500, 0.5500, 0.5900, 0.6000, 0.6100, 0.6300, 0.6400, 0.6400, 0.6700, 0.6800, 0.6800,
394 0.6900, 0.7100, 0.7100, 0.7100, 0.7200, 0.7300, 0.7400, 0.7400, 0.7500, 0.7500, 0.7600, 0.7600, 0.7700, 0.7700,
395 0.7800, 0.7800, 0.7800, 0.7900, 0.7900, 0.7900, 0.8000, 0.8100, 0.8100, 0.8100, 0.8200, 0.8200, 0.8300, 0.8300,
396 0.8400, 0.8400, 0.8400, 0.8400, 0.8400, 0.8400, 0.8500, 0.8600, 0.0550, 0.1200, 0.1400, 0.1500, 0.1700, 0.1800,
397 0.2300, 0.2700, 0.3100, 0.3600, 0.3800, 0.4100, 0.4200, 0.4600, 0.4800, 0.5000, 0.5100, 0.5400, 0.5500, 0.5800,
398 0.5900, 0.6000, 0.6200, 0.6400, 0.6400, 0.6600, 0.6700, 0.6800, 0.6900, 0.7000, 0.7000, 0.7000, 0.7200, 0.7300,
399 0.7300, 0.7400, 0.7500, 0.7500, 0.7500, 0.7600, 0.7700, 0.7700, 0.7700, 0.7800, 0.7800, 0.7800, 0.7900, 0.7900,
400 0.8000, 0.8000, 0.8100, 0.8100, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400, 0.8400, 0.8400,
401 0.8500, 0.8600, 0.0600, 0.1200, 0.1300, 0.1500, 0.1600, 0.1700, 0.2100, 0.2600, 0.3000, 0.3600, 0.3700, 0.4000,
402 0.4200, 0.4600, 0.4700, 0.5000, 0.5100, 0.5400, 0.5400, 0.5800, 0.5900, 0.6000, 0.6200, 0.6300, 0.6300, 0.6600,
403 0.6700, 0.6700, 0.6900, 0.7000, 0.7000, 0.7000, 0.7100, 0.7200, 0.7300, 0.7400, 0.7400, 0.7500, 0.7500, 0.7600,
404 0.7600, 0.7700, 0.7700, 0.7800, 0.7800, 0.7800, 0.7900, 0.7900, 0.8000, 0.8000, 0.8100, 0.8100, 0.8100, 0.8200,
405 0.8200, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400, 0.8400, 0.8400, 0.8500, 0.8600, 0.0650, 0.1100, 0.1200, 0.1400,
406 0.1500, 0.1600, 0.2000, 0.2600, 0.3000, 0.3500, 0.3600, 0.4000, 0.4100, 0.4500, 0.4600, 0.5000, 0.5100, 0.5300,
407 0.5400, 0.5800, 0.5900, 0.5900, 0.6100, 0.6300, 0.6300, 0.6600, 0.6700, 0.6700, 0.6800, 0.7000, 0.7000, 0.7000,
408 0.7100, 0.7200, 0.7300, 0.7300, 0.7400, 0.7500, 0.7500, 0.7600, 0.7600, 0.7700, 0.7700, 0.7700, 0.7800, 0.7800,
409 0.7900, 0.7900, 0.8000, 0.8000, 0.8100, 0.8100, 0.8100, 0.8200, 0.8200, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400,
410 0.8400, 0.8400, 0.8500, 0.8500, 0.0700, 0.1100, 0.1200, 0.1300, 0.1400, 0.1500, 0.2000, 0.2500, 0.3000, 0.3500,
411 0.3600, 0.3900, 0.4100, 0.4500, 0.4600, 0.4900, 0.5000, 0.5300, 0.5400, 0.5800, 0.5800, 0.5900, 0.6100, 0.6300,
412 0.6300, 0.6600, 0.6600, 0.6700, 0.6800, 0.6900, 0.6900, 0.7000, 0.7100, 0.7200, 0.7300, 0.7300, 0.7400, 0.7400,
413 0.7500, 0.7600, 0.7600, 0.7700, 0.7700, 0.7700, 0.7800, 0.7800, 0.7800, 0.7900, 0.8000, 0.8000, 0.8100, 0.8100,
414 0.8100, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400, 0.8400, 0.8500, 0.8500, 0.0750, 0.1000,
415 0.1100, 0.1200, 0.1300, 0.1500, 0.1900, 0.2400, 0.2900, 0.3400, 0.3600, 0.3900, 0.4100, 0.4400, 0.4500, 0.4900,
416 0.5000, 0.5200, 0.5400, 0.5700, 0.5800, 0.5900, 0.6100, 0.6300, 0.6300, 0.6500, 0.6600, 0.6700, 0.6800, 0.6900,
417 0.6900, 0.6900, 0.7100, 0.7200, 0.7300, 0.7300, 0.7400, 0.7400, 0.7500, 0.7600, 0.7600, 0.7700, 0.7700, 0.7700,
418 0.7800, 0.7800, 0.7800, 0.7900, 0.8000, 0.8000, 0.8000, 0.8100, 0.8100, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300,
419 0.8400, 0.8400, 0.8400, 0.8400, 0.8400, 0.8500, 0.0800, 0.1000, 0.1000, 0.1200, 0.1300, 0.1400, 0.1900, 0.2300,
420 0.2900, 0.3400, 0.3500, 0.3800, 0.4000, 0.4400, 0.4500, 0.4900, 0.5000, 0.5200, 0.5300, 0.5700, 0.5800, 0.5900,
421 0.6100, 0.6200, 0.6200, 0.6500, 0.6600, 0.6700, 0.6800, 0.6900, 0.6900, 0.6900, 0.7100, 0.7200, 0.7300, 0.7300,
422 0.7400, 0.7400, 0.7500, 0.7600, 0.7600, 0.7600, 0.7700, 0.7700, 0.7700, 0.7800, 0.7800, 0.7900, 0.7900, 0.8000,
423 0.8000, 0.8000, 0.8100, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300, 0.8400, 0.8400, 0.8400, 0.8400, 0.8400, 0.8500,
424 0.0850, 0.1000, 0.1000, 0.1100, 0.1200, 0.1400, 0.1800, 0.2300, 0.2800, 0.3300, 0.3500, 0.3800, 0.4000, 0.4400,
425 0.4500, 0.4900, 0.5000, 0.5100, 0.5300, 0.5700, 0.5800, 0.5900, 0.6100, 0.6200, 0.6200, 0.6500, 0.6600, 0.6700,
426 0.6800, 0.6800, 0.6900, 0.6900, 0.7100, 0.7200, 0.7300, 0.7300, 0.7400, 0.7400, 0.7500, 0.7600, 0.7600, 0.7600,
427 0.7600, 0.7700, 0.7700, 0.7800, 0.7800, 0.7900, 0.7900, 0.8000, 0.8000, 0.8000, 0.8100, 0.8200, 0.8200, 0.8200,
428 0.8300, 0.8300, 0.8300, 0.8300, 0.8300, 0.8400, 0.8400, 0.8500, 0.0900, 0.0900, 0.0900, 0.1100, 0.1200, 0.1300,
429 0.1800, 0.2200, 0.2800, 0.3300, 0.3500, 0.3800, 0.4000, 0.4300, 0.4500, 0.4900, 0.5000, 0.5100, 0.5300, 0.5600,
430 0.5800, 0.5900, 0.6100, 0.6200, 0.6200, 0.6500, 0.6600, 0.6600, 0.6700, 0.6700, 0.6900, 0.6900, 0.7100, 0.7100,
431 0.7300, 0.7300, 0.7300, 0.7400, 0.7500, 0.7600, 0.7600, 0.7600, 0.7600, 0.7700, 0.7700, 0.7800, 0.7800, 0.7900,
432 0.7900, 0.8000, 0.8000, 0.8000, 0.8100, 0.8100, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300, 0.8300, 0.8300, 0.8400,
433 0.8400, 0.8500, 0.0950, 0.0900, 0.0900, 0.1000, 0.1100, 0.1300, 0.1800, 0.2200, 0.2800, 0.3300, 0.3400, 0.3800,
434 0.4000, 0.4300, 0.4500, 0.4800, 0.4900, 0.5100, 0.5300, 0.5600, 0.5700, 0.5900, 0.6100, 0.6200, 0.6200, 0.6500,
435 0.6500, 0.6600, 0.6700, 0.6700, 0.6800, 0.6900, 0.7100, 0.7100, 0.7200, 0.7300, 0.7300, 0.7400, 0.7500, 0.7500,
436 0.7600, 0.7600, 0.7600, 0.7700, 0.7700, 0.7800, 0.7800, 0.7900, 0.7900, 0.8000, 0.8000, 0.8000, 0.8100, 0.8100,
437 0.8200, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300, 0.8300, 0.8400, 0.8400, 0.8500, 0.1000, 0.0900, 0.0900, 0.1000,
438 0.1100, 0.1200, 0.1700, 0.2200, 0.2700, 0.3200, 0.3400, 0.3700, 0.4000, 0.4300, 0.4400, 0.4800, 0.4900, 0.5100,
439 0.5300, 0.5600, 0.5700, 0.5800, 0.6100, 0.6200, 0.6200, 0.6400, 0.6500, 0.6600, 0.6700, 0.6700, 0.6800, 0.6900,
440 0.7100, 0.7100, 0.7200, 0.7300, 0.7300, 0.7400, 0.7500, 0.7500, 0.7600, 0.7600, 0.7600, 0.7700, 0.7700, 0.7800,
441 0.7800, 0.7900, 0.7900, 0.8000, 0.8000, 0.8000, 0.8100, 0.8100, 0.8200, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300,
442 0.8300, 0.8300, 0.8400, 0.8500, 0.1200, 0.0800, 0.0800, 0.0900, 0.0900, 0.1100, 0.1600, 0.2100, 0.2600, 0.3100,
443 0.3300, 0.3700, 0.3900, 0.4200, 0.4400, 0.4800, 0.4900, 0.5100, 0.5300, 0.5500, 0.5700, 0.5800, 0.6000, 0.6100,
444 0.6100, 0.6400, 0.6400, 0.6600, 0.6700, 0.6700, 0.6800, 0.6800, 0.7000, 0.7100, 0.7200, 0.7200, 0.7300, 0.7400,
445 0.7400, 0.7500, 0.7600, 0.7600, 0.7600, 0.7700, 0.7700, 0.7800, 0.7800, 0.7800, 0.7900, 0.7900, 0.8000, 0.8000,
446 0.8100, 0.8100, 0.8100, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300, 0.8300, 0.8300, 0.8400, 0.8400, 0.1500, 0.0700,
447 0.0700, 0.0800, 0.0800, 0.1000, 0.1600, 0.2000, 0.2500, 0.2900, 0.3200, 0.3500, 0.3900, 0.4100, 0.4400, 0.4700,
448 0.4900, 0.5000, 0.5200, 0.5500, 0.5600, 0.5700, 0.6000, 0.6100, 0.6100, 0.6400, 0.6400, 0.6600, 0.6600, 0.6600,
449 0.6800, 0.6800, 0.7000, 0.7100, 0.7100, 0.7200, 0.7300, 0.7300, 0.7400, 0.7500, 0.7600, 0.7600, 0.7600, 0.7600,
450 0.7700, 0.7700, 0.7800, 0.7800, 0.7900, 0.7900, 0.7900, 0.8000, 0.8100, 0.8100, 0.8100, 0.8100, 0.8200, 0.8200,
451 0.8300, 0.8300, 0.8300, 0.8300, 0.8400, 0.8400, 0.2000, 0.0500, 0.0500, 0.0600, 0.0600, 0.1000, 0.1500, 0.1900,
452 0.2400, 0.2700, 0.3100, 0.3400, 0.3800, 0.4100, 0.4300, 0.4600, 0.4800, 0.5000, 0.5100, 0.5400, 0.5500, 0.5700,
453 0.5900, 0.6000, 0.6100, 0.6400, 0.6400, 0.6500, 0.6600, 0.6600, 0.6800, 0.6800, 0.7000, 0.7100, 0.7100, 0.7100,
454 0.7300, 0.7300, 0.7400, 0.7500, 0.7600, 0.7600, 0.7600, 0.7600, 0.7600, 0.7700, 0.7700, 0.7800, 0.7800, 0.7900,
455 0.7900, 0.8000, 0.8000, 0.8100, 0.8100, 0.8100, 0.8100, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300, 0.8300, 0.8400,
456 0.2500, 0.0400, 0.0400, 0.0500, 0.0600, 0.0900, 0.1400, 0.1800, 0.2300, 0.2500, 0.3000, 0.3400, 0.3700, 0.4100,
457 0.4300, 0.4600, 0.4700, 0.5000, 0.5100, 0.5400, 0.5500, 0.5700, 0.5900, 0.6000, 0.6000, 0.6300, 0.6400, 0.6500,
458 0.6600, 0.6600, 0.6800, 0.6800, 0.7000, 0.7000, 0.7100, 0.7100, 0.7300, 0.7300, 0.7400, 0.7500, 0.7500, 0.7600,
459 0.7600, 0.7600, 0.7600, 0.7700, 0.7700, 0.7700, 0.7800, 0.7800, 0.7900, 0.8000, 0.8000, 0.8100, 0.8100, 0.8100,
460 0.8100, 0.8200, 0.8200, 0.8200, 0.8200, 0.8300, 0.8300, 0.8400, 0.2600, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
461 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
462 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
463 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
464 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
465 0.0000, 0.0000, 1.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
466 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
467 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
468 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
469 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
470};
471static const float inv_ratio_tbl_10[] = {
472 0.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
473 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
474 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.0100,
475 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
476 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
477 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.0110, 0.9635,
478 0.9662, 0.9690, 0.9714, 0.9732, 0.9781, 0.9777, 0.9786, 0.9801, 0.9816, 0.9761, 0.9706, 0.9762, 0.9817, 0.9811,
479 0.9842, 0.9842, 0.9842, 0.9868, 0.9894, 0.9874, 0.9876, 0.9869, 0.9862, 0.9891, 0.9919, 0.9880, 0.9922, 0.9902,
480 0.9882, 0.9920, 0.9958, 0.9958, 0.9958, 0.9958, 0.9958, 0.9958, 0.9958, 0.9958, 0.9958, 0.0160, 0.7998, 0.8149,
481 0.8300, 0.8406, 0.8515, 0.8781, 0.8800, 0.8837, 0.8905, 0.8973, 0.8969, 0.8965, 0.9054, 0.9142, 0.9204, 0.9141,
482 0.9221, 0.9301, 0.9342, 0.9382, 0.9378, 0.9415, 0.9422, 0.9429, 0.9476, 0.9523, 0.9463, 0.9591, 0.9553, 0.9516,
483 0.9635, 0.9755, 0.9755, 0.9755, 0.9755, 0.9755, 0.9755, 0.9755, 0.9755, 0.9755, 0.0210, 0.6824, 0.7028, 0.7232,
484 0.7435, 0.7621, 0.7951, 0.7999, 0.8191, 0.8289, 0.8387, 0.8424, 0.8461, 0.8566, 0.8670, 0.8779, 0.8682, 0.8793,
485 0.8904, 0.8949, 0.8995, 0.9033, 0.9089, 0.9109, 0.9129, 0.9161, 0.9194, 0.9196, 0.9296, 0.9277, 0.9257, 0.9399,
486 0.9541, 0.9541, 0.9541, 0.9541, 0.9541, 0.9541, 0.9541, 0.9541, 0.9541, 0.0260, 0.5879, 0.6139, 0.6400, 0.6721,
487 0.7004, 0.7297, 0.7405, 0.7750, 0.7873, 0.7996, 0.8024, 0.8053, 0.8202, 0.8352, 0.8446, 0.8439, 0.8546, 0.8653,
488 0.8704, 0.8754, 0.8811, 0.8841, 0.8879, 0.8917, 0.8953, 0.8988, 0.8982, 0.9055, 0.9062, 0.9069, 0.9197, 0.9325,
489 0.9325, 0.9325, 0.9325, 0.9325, 0.9325, 0.9325, 0.9325, 0.9325, 0.0310, 0.5198, 0.5514, 0.5830, 0.6154, 0.6500,
490 0.6794, 0.7006, 0.7422, 0.7556, 0.7690, 0.7723, 0.7756, 0.7928, 0.8099, 0.8188, 0.8263, 0.8351, 0.8439, 0.8494,
491 0.8549, 0.8622, 0.8667, 0.8711, 0.8756, 0.8790, 0.8824, 0.8821, 0.8890, 0.8910, 0.8929, 0.9026, 0.9123, 0.9123,
492 0.9123, 0.9123, 0.9123, 0.9123, 0.9123, 0.9123, 0.9123, 0.0360, 0.4785, 0.5086, 0.5388, 0.5707, 0.5993, 0.6428,
493 0.6693, 0.7145, 0.7290, 0.7435, 0.7490, 0.7546, 0.7721, 0.7896, 0.7973, 0.8094, 0.8176, 0.8258, 0.8316, 0.8374,
494 0.8459, 0.8527, 0.8577, 0.8627, 0.8653, 0.8680, 0.8696, 0.8767, 0.8796, 0.8824, 0.8886, 0.8947, 0.8947, 0.8947,
495 0.8947, 0.8947, 0.8947, 0.8947, 0.8947, 0.8947, 0.0410, 0.4392, 0.4715, 0.5039, 0.5337, 0.5535, 0.6109, 0.6434,
496 0.6894, 0.7046, 0.7198, 0.7277, 0.7356, 0.7542, 0.7729, 0.7787, 0.7941, 0.8023, 0.8104, 0.8165, 0.8226, 0.8319,
497 0.8408, 0.8462, 0.8516, 0.8535, 0.8555, 0.8588, 0.8670, 0.8706, 0.8743, 0.8779, 0.8814, 0.8814, 0.8814, 0.8814,
498 0.8814, 0.8814, 0.8814, 0.8814, 0.8814, 0.0460, 0.4110, 0.4431, 0.4751, 0.5024, 0.5223, 0.5831, 0.6201, 0.6669,
499 0.6821, 0.6973, 0.7075, 0.7178, 0.7378, 0.7578, 0.7634, 0.7802, 0.7888, 0.7975, 0.8039, 0.8103, 0.8198, 0.8309,
500 0.8365, 0.8420, 0.8434, 0.8448, 0.8495, 0.8597, 0.8632, 0.8667, 0.8700, 0.8732, 0.8732, 0.8732, 0.8732, 0.8732,
501 0.8732, 0.8732, 0.8732, 0.8732, 0.0510, 0.3880, 0.4187, 0.4494, 0.4754, 0.4960, 0.5589, 0.5991, 0.6470, 0.6617,
502 0.6765, 0.6895, 0.7024, 0.7233, 0.7442, 0.7507, 0.7672, 0.7768, 0.7865, 0.7933, 0.8001, 0.8094, 0.8229, 0.8283,
503 0.8338, 0.8348, 0.8359, 0.8412, 0.8537, 0.8567, 0.8597, 0.8632, 0.8666, 0.8666, 0.8666, 0.8666, 0.8666, 0.8666,
504 0.8666, 0.8666, 0.8666, 0.0560, 0.3673, 0.3969, 0.4265, 0.4515, 0.4734, 0.5380, 0.5803, 0.6297, 0.6440, 0.6582,
505 0.6736, 0.6890, 0.7104, 0.7318, 0.7394, 0.7561, 0.7665, 0.7769, 0.7842, 0.7916, 0.8001, 0.8166, 0.8216, 0.8266,
506 0.8276, 0.8286, 0.8337, 0.8480, 0.8506, 0.8532, 0.8570, 0.8608, 0.8608, 0.8608, 0.8608, 0.8608, 0.8608, 0.8608,
507 0.8608, 0.8608, 0.0610, 0.3492, 0.3777, 0.4062, 0.4303, 0.4535, 0.5199, 0.5639, 0.6151, 0.6290, 0.6429, 0.6598,
508 0.6768, 0.6989, 0.7210, 0.7294, 0.7469, 0.7578, 0.7686, 0.7763, 0.7840, 0.7918, 0.8109, 0.8154, 0.8199, 0.8214,
509 0.8229, 0.8267, 0.8426, 0.8449, 0.8472, 0.8514, 0.8556, 0.8556, 0.8556, 0.8556, 0.8556, 0.8556, 0.8556, 0.8556,
510 0.8556, 0.0660, 0.3337, 0.3611, 0.3885, 0.4118, 0.4352, 0.5042, 0.5498, 0.6030, 0.6170, 0.6311, 0.6483, 0.6655,
511 0.6885, 0.7115, 0.7207, 0.7380, 0.7496, 0.7612, 0.7693, 0.7774, 0.7845, 0.8055, 0.8095, 0.8136, 0.8159, 0.8181,
512 0.8208, 0.8375, 0.8396, 0.8417, 0.8464, 0.8510, 0.8510, 0.8510, 0.8510, 0.8510, 0.8510, 0.8510, 0.8510, 0.8510,
513 0.0710, 0.3209, 0.3470, 0.3731, 0.3962, 0.4182, 0.4902, 0.5370, 0.5923, 0.6072, 0.6221, 0.6386, 0.6550, 0.6792,
514 0.7035, 0.7131, 0.7294, 0.7419, 0.7545, 0.7630, 0.7715, 0.7782, 0.8003, 0.8041, 0.8078, 0.8107, 0.8137, 0.8160,
515 0.8327, 0.8347, 0.8367, 0.8419, 0.8470, 0.8470, 0.8470, 0.8470, 0.8470, 0.8470, 0.8470, 0.8470, 0.8470, 0.0760,
516 0.3100, 0.3350, 0.3601, 0.3834, 0.4028, 0.4777, 0.5248, 0.5828, 0.5983, 0.6138, 0.6296, 0.6454, 0.6709, 0.6965,
517 0.7059, 0.7213, 0.7348, 0.7483, 0.7574, 0.7664, 0.7727, 0.7955, 0.7989, 0.8024, 0.8060, 0.8097, 0.8120, 0.8283,
518 0.8303, 0.8323, 0.8379, 0.8436, 0.8436, 0.8436, 0.8436, 0.8436, 0.8436, 0.8436, 0.8436, 0.8436, 0.0810, 0.3003,
519 0.3245, 0.3488, 0.3723, 0.3894, 0.4660, 0.5133, 0.5743, 0.5901, 0.6060, 0.6214, 0.6368, 0.6633, 0.6899, 0.6992,
520 0.7138, 0.7281, 0.7424, 0.7522, 0.7620, 0.7681, 0.7909, 0.7941, 0.7974, 0.8017, 0.8060, 0.8084, 0.8243, 0.8263,
521 0.8284, 0.8345, 0.8406, 0.8406, 0.8406, 0.8406, 0.8406, 0.8406, 0.8406, 0.8406, 0.8406, 0.0860, 0.2912, 0.3149,
522 0.3386, 0.3621, 0.3783, 0.4547, 0.5026, 0.5664, 0.5826, 0.5987, 0.6139, 0.6291, 0.6564, 0.6837, 0.6929, 0.7069,
523 0.7218, 0.7367, 0.7474, 0.7582, 0.7637, 0.7865, 0.7898, 0.7930, 0.7978, 0.8027, 0.8052, 0.8206, 0.8228, 0.8250,
524 0.8315, 0.8381, 0.8381, 0.8381, 0.8381, 0.8381, 0.8381, 0.8381, 0.8381, 0.8381, 0.0910, 0.2825, 0.3055, 0.3285,
525 0.3528, 0.3685, 0.4439, 0.4928, 0.5591, 0.5755, 0.5920, 0.6070, 0.6220, 0.6499, 0.6778, 0.6871, 0.7009, 0.7160,
526 0.7312, 0.7430, 0.7548, 0.7596, 0.7824, 0.7858, 0.7891, 0.7944, 0.7997, 0.8023, 0.8172, 0.8196, 0.8220, 0.8289,
527 0.8359, 0.8359, 0.8359, 0.8359, 0.8359, 0.8359, 0.8359, 0.8359, 0.8359, 0.0960, 0.2736, 0.2958, 0.3179, 0.3441,
528 0.3595, 0.4337, 0.4839, 0.5519, 0.5688, 0.5856, 0.6008, 0.6160, 0.6442, 0.6723, 0.6818, 0.6957, 0.7109, 0.7262,
529 0.7388, 0.7515, 0.7558, 0.7786, 0.7821, 0.7857, 0.7913, 0.7969, 0.7998, 0.8142, 0.8167, 0.8192, 0.8265, 0.8339,
530 0.8339, 0.8339, 0.8339, 0.8339, 0.8339, 0.8339, 0.8339, 0.8339, 0.1200, 0.2294, 0.2497, 0.2700, 0.2989, 0.3241,
531 0.3954, 0.4541, 0.5183, 0.5395, 0.5607, 0.5801, 0.5995, 0.6252, 0.6509, 0.6622, 0.6812, 0.6958, 0.7104, 0.7239,
532 0.7373, 0.7418, 0.7635, 0.7684, 0.7733, 0.7802, 0.7872, 0.7907, 0.8051, 0.8069, 0.8087, 0.8171, 0.8254, 0.8254,
533 0.8254, 0.8254, 0.8254, 0.8254, 0.8254, 0.8254, 0.8254, 0.1500, 0.1736, 0.2027, 0.2318, 0.2595, 0.2893, 0.3658,
534 0.4295, 0.4860, 0.5119, 0.5378, 0.5621, 0.5863, 0.6100, 0.6338, 0.6497, 0.6696, 0.6850, 0.7004, 0.7120, 0.7236,
535 0.7327, 0.7517, 0.7576, 0.7634, 0.7718, 0.7802, 0.7850, 0.7985, 0.8008, 0.8032, 0.8103, 0.8175, 0.8175, 0.8175,
536 0.8175, 0.8175, 0.8175, 0.8175, 0.8175, 0.8175, 0.2000, 0.0986, 0.1392, 0.1797, 0.2181, 0.2581, 0.3283, 0.4013,
537 0.4499, 0.4799, 0.5099, 0.5390, 0.5680, 0.5921, 0.6162, 0.6369, 0.6574, 0.6725, 0.6876, 0.6995, 0.7114, 0.7234,
538 0.7417, 0.7475, 0.7534, 0.7631, 0.7728, 0.7797, 0.7904, 0.7943, 0.7983, 0.8034, 0.8085, 0.8085, 0.8085, 0.8085,
539 0.8085, 0.8085, 0.8085, 0.8085, 0.8085, 0.2500, 0.0394, 0.0909, 0.1424, 0.1873, 0.2355, 0.2974, 0.3791, 0.4209,
540 0.4555, 0.4900, 0.5217, 0.5535, 0.5799, 0.6062, 0.6283, 0.6476, 0.6621, 0.6766, 0.6909, 0.7051, 0.7173, 0.7305,
541 0.7389, 0.7474, 0.7578, 0.7681, 0.7749, 0.7851, 0.7894, 0.7937, 0.7975, 0.8013, 0.8013, 0.8013, 0.8013, 0.8013,
542 0.8013, 0.8013, 0.8013, 0.8013, 0.3000, 0.0210, 0.0711, 0.1212, 0.1689, 0.2136, 0.2731, 0.3608, 0.3978, 0.4363,
543 0.4748, 0.5090, 0.5432, 0.5711, 0.5989, 0.6215, 0.6396, 0.6545, 0.6693, 0.6852, 0.7011, 0.7122, 0.7202, 0.7315,
544 0.7428, 0.7534, 0.7641, 0.7707, 0.7816, 0.7858, 0.7899, 0.7927, 0.7954, 0.7954, 0.7954, 0.7954, 0.7954, 0.7954,
545 0.7954, 0.7954, 0.7954, 0.3500, 0.0113, 0.0600, 0.1086, 0.1539, 0.1971, 0.2541, 0.3459, 0.3831, 0.4243, 0.4654,
546 0.5010, 0.5367, 0.5648, 0.5930, 0.6164, 0.6340, 0.6490, 0.6639, 0.6814, 0.6988, 0.7084, 0.7147, 0.7270, 0.7393,
547 0.7502, 0.7610, 0.7668, 0.7787, 0.7830, 0.7873, 0.7899, 0.7925, 0.7925, 0.7925, 0.7925, 0.7925, 0.7925, 0.7925,
548 0.7925, 0.7925, 0.4000, 0.0112, 0.0555, 0.0997, 0.1425, 0.1848, 0.2390, 0.3341, 0.3721, 0.4153, 0.4586, 0.4950,
549 0.5315, 0.5593, 0.5872, 0.6129, 0.6314, 0.6454, 0.6593, 0.6782, 0.6971, 0.7068, 0.7139, 0.7253, 0.7366, 0.7479,
550 0.7593, 0.7657, 0.7765, 0.7810, 0.7855, 0.7889, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924,
551 0.7924, 0.4500, 0.0111, 0.0520, 0.0928, 0.1337, 0.1751, 0.2291, 0.3245, 0.3613, 0.4086, 0.4558, 0.4914, 0.5270,
552 0.5544, 0.5818, 0.6101, 0.6298, 0.6430, 0.6562, 0.6760, 0.6958, 0.7057, 0.7134, 0.7240, 0.7345, 0.7464, 0.7584,
553 0.7651, 0.7758, 0.7799, 0.7841, 0.7882, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924,
554 0.4990, 0.0086, 0.0484, 0.0881, 0.1279, 0.1677, 0.2232, 0.3162, 0.3570, 0.4060, 0.4550, 0.4888, 0.5227, 0.5499,
555 0.5771, 0.6086, 0.6285, 0.6417, 0.6550, 0.6752, 0.6954, 0.7053, 0.7132, 0.7232, 0.7331, 0.7455, 0.7579, 0.7649,
556 0.7758, 0.7796, 0.7834, 0.7879, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.7924, 0.5000,
557 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
558 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
559 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, 0.0000,
560 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
561 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
562 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
563};
564/* cdf of normal distribution -----------------------------------------------*/
565static double cdf_norm(double x) {
566 static const double cdf_tbl[] = {
567 /* cdf(Normal(0,1),0:0.1:7) */
568 0.500000000000, 0.539827837277, 0.579259709439, 0.617911422189, 0.655421741610, 0.691462461274, 0.725746882250,
569 0.758036347777, 0.788144601417, 0.815939874653, 0.841344746069, 0.864333939054, 0.884930329778, 0.903199515414,
570 0.919243340766, 0.933192798731, 0.945200708300, 0.955434537241, 0.964069680887, 0.971283440184, 0.977249868052,
571 0.982135579437, 0.986096552487, 0.989275889978, 0.991802464075, 0.993790334674, 0.995338811976, 0.996533026197,
572 0.997444869670, 0.998134186700, 0.998650101968, 0.999032396787, 0.999312862062, 0.999516575858, 0.999663070734,
573 0.999767370921, 0.999840891410, 0.999892200267, 0.999927651956, 0.999951903656, 0.999968328758, 0.999979342493,
574 0.999986654251, 0.999991460095, 0.999994587456, 0.999996602327, 0.999997887545, 0.999998699193, 0.999999206672,
575 0.999999520817, 0.999999713348, 0.999999830173, 0.999999900356, 0.999999942099, 0.999999966680, 0.999999981010,
576 0.999999989282, 0.999999994010, 0.999999996684, 0.999999998182, 0.999999999013, 0.999999999470, 0.999999999718,
577 0.999999999851, 0.999999999922, 0.999999999960, 0.999999999979, 0.999999999990, 0.999999999995, 0.999999999997,
578 0.999999999999};
579 double a = x / 0.1;
580 int i = (int)floor(a);
581
582 if (i < 0) {
583 return 0.0; /* x < 0.0 */
584 }
585 if (i > 69) {
586 return 1.0;
587 }
588 return (1.0 - a + i) * cdf_tbl[i] + (a - i) * cdf_tbl[i + 1];
589}
LD
static int LD(int n, const double *Q, double *L, double *D)
Definition mrtk_lambda.c:40
lambda
int lambda(int n, int m, const double *a, const double *Q, double *F, double *s)
Integer least-square estimation (LAMBDA/MLAMBDA).
Definition mrtk_lambda.c:225
lambda_search
int lambda_search(int n, int m, const double *a, const double *Q, double *F, double *s)
MLAMBDA search for integer least-squares.
Definition mrtk_lambda.c:302
SWAP
#define SWAP(x, y)
Definition mrtk_lambda.c:31
LOOPMAX
#define LOOPMAX
Definition mrtk_lambda.c:27
SGN
#define SGN(x)
Definition mrtk_lambda.c:29
gauss
static void gauss(int n, double *L, double *Z, int i, int j)
Definition mrtk_lambda.c:70
inv_ratio_tbl_1
static const float inv_ratio_tbl_1[]
Definition mrtk_lambda.c:327
ROUND
#define ROUND(x)
Definition mrtk_lambda.c:30
inv_ratio_tbl_10
static const float inv_ratio_tbl_10[]
Definition mrtk_lambda.c:471
search
static int search(int n, int m, const double *L, const double *D, const double *zs, double *zn, double *s)
Definition mrtk_lambda.c:129
perm
static void perm(int n, double *L, double *D, int j, double del, double *Z)
Definition mrtk_lambda.c:83
cdf_norm
static double cdf_norm(double x)
Definition mrtk_lambda.c:565
reduction
static void reduction(int n, double *L, double *D, double *Z)
Definition mrtk_lambda.c:106
lambda_reduction
int lambda_reduction(int n, const double *Q, double *Z)
LAMBDA reduction for integer least-squares.
Definition mrtk_lambda.c:263
mrtk_lambda.h
MRTKLIB Integer Ambiguity Resolution Module — LAMBDA method.
mat
double * mat(int n, int m)
Allocate memory for a double matrix.
Definition mrtk_mat.c:76
zeros
double * zeros(int n, int m)
Allocate a new zero matrix.
Definition mrtk_mat.c:108
solve
int solve(const char *tr, const double *A, const double *Y, int n, int m, double *X)
Solve linear equation (X = A\Y or X = A'\Y).
Definition mrtk_mat.c:443
matmul
void matmul(const char *tr, int n, int k, int m, double alpha, const double *A, const double *B, double beta, double *C)
Multiply matrix by matrix (C = alpha*A*B + beta*C).
Definition mrtk_mat.c:310
eye
double * eye(int n)
Generate a new identity matrix.
Definition mrtk_mat.c:129
mrtk_mat.h
MRTKLIB Matrix Module — Matrix/vector operations and estimation.
free
#define free(x)
Definition toml.c:53