/src/reduction.c
C | 89 lines | 68 code | 17 blank | 4 comment | 10 complexity | 1c36dc3b582d12483f91e111318a1750 MD5 | raw file
1#include "reduction.h"
2
3float** bezier_reduction_matrix(int n, int m)
4{
5 assert(n > m);
6
7 float (*psi)[n+1] = malloc(sizeof(float) * (m + 1) * (n + 1));
8
9 float C = factrl(n)*up_power(1-n, m)
10 / (factrl(m)*up_power(m+2, n));
11
12 // 0,0
13 {
14 psi[0][0] = 0;
15 for(int h = 0; h <= m; ++h)
16 {
17 float t = bico(m,h) * up_power(m+2,h) / ((h-n)*factrl(h));
18 if(h % 2 != 0)
19 t *= -1;
20 psi[0][0] += t;
21 }
22 psi[0][0] *= -n*C;
23 }
24
25 // 0,1
26 {
27 psi[0][1] = (m+1)*C;
28 if(m%2 != 0)
29 psi[0][1] *= -1;
30 }
31
32 // first row
33 {
34 for(int j = 1; j <= n - 1; ++j)
35 {
36 float F = ((float)((j-1)*(n-j+1)))/((float) ((n-j)*(j+1)));
37 float E = 1. + F - ((float)(m*(m+2)))/((float) ((n-j)*(j+1)));
38
39 psi[0][j+1] = E*psi[0][j] - F*psi[0][j-1];
40 }
41 }
42
43 // rest of matrix
44 {
45 for(int i = 0; i <= m - 1; ++i)
46 {
47 for(int j = 0; j <= n; ++j)
48 {
49 float Anj = j*(j-n-1);
50 float Ami = i*(i-m-1);
51 float Bnj = (j-n)*(j+1);
52 float Bmi = (i-m)*(i+1);
53 float Cnj = Anj + Bnj;
54 float Cmi = Ami + Bmi;
55
56 float pp0m1 = ((j - 1) < 0) ? 0 : psi[i][j-1];
57 float pp0p0 = psi[i][j];
58 float pm1p0 = ((i - 1) < 0) ? 0 : psi[i-1][j];
59 float pp0p1 = ((j + 1) > n) ? 0 : psi[i][j+1];
60
61 psi[i+1][j] = (Anj*pp0m1 + (Cmi-Cnj)*pp0p0 +
62 Bnj*pp0p1 - Ami*pm1p0) / Bmi;
63 }
64 }
65 }
66
67 return (float**)psi;
68}
69
70Bezier* bezier_degree_reduction(Bezier* c, int m, float **reduction_matrix)
71{
72 int n = c->n;
73
74 float (*psi)[n+1] = (float(*)[])reduction_matrix;
75 float (*e) = malloc(sizeof(float) * (n + 1));
76
77 for(int i = 0; i <= m; ++i)
78 {
79 e[i] = 0;
80 for(int j = 0; j <= n; ++j)
81 e[i] += c->c[j]*psi[i][j];
82 }
83
84 Bezier* low = bezier_create_with_coeffs(m, e);
85 low->dom = interval_copy(c->dom);
86
87 return low;
88}
89