[prev in list] [next in list] [prev in thread] [next in thread]
List: kde-commits
Subject: koffice/krita/plugins/viewplugins/painterlymixer
From: Emanuele Tamponi <emanuele () valinor ! it>
Date: 2007-08-07 15:55:41
Message-ID: 1186502141.097148.15448.nullmailer () svn ! kde ! org
[Download RAW message or body]
SVN commit 697377 by tamponi:
Use GMP++ everywhere.
M +30 -30 IlluminantD65.ill
M +1 -1 ipg/compile
M +130 -160 ipg/ipg.cpp
--- trunk/koffice/krita/plugins/viewplugins/painterlymixer/IlluminantD65.ill \
#697376:697377 @@ -1,30 +1,30 @@
-1.871929361865446091286087571537431230306314233757802867330610752e-06
-0.0002360339403767762923554329637220589033219653174455743283033370972
-0.003797070661949299869431499352898096155684015684528276324272155762
-0.05017762889547899757188573546806686920263018691912293434143066406
-0.1098932324552358724643603797577195990697873639874160289764404297
-0.00764154200729675839760406163073125362927839887561276555061340332
-0.01601480499693024116931398120877716451104788575321435928344726562
-0.1915942900317865673785657637528956342976016458123922348022460938
-0.1779143309921572162413854270579527394602337153628468513488769531
-0.3810041601574760195328666484426705096666410099714994430541992188
-7.5029084263711804391353376235471032192764218038405488186981529e-07
-9.22123262305353224317144057188815747805321620944596361368894577e-05
-0.0003975693610413827529699514157050799223469539356301538646221160889
-0.005504467680665073033517540766410003350017632328672334551811218262
-0.02930374594476542653128444283067155495814404275733977556228637695
-0.1998636693728007917929788417321645965785137377679347991943359375
-0.00624253998238867203085696033770801705031772144138813018798828125
-0.4167196173230689423578210900389606763383198995143175125122070312
-0.07706609115986708507655420299919946103273105109110474586486816406
-0.2648093365583294538621612207451860854234837461262941360473632812
-9.999999999999999653057388335469287129029766072807834803722132479e-65
-9.999999999999999653057388335469287129029766072807834803722132479e-65
-0.01711914340041561673426537184411477809931056981440633535385131836
-0.2421760160212150506982156997204391757350094849243760108947753906
-0.598189061596961770494139398746114011373720131814479827880859375
-0.2240890811920233766888366916392172356609080452471971511840820312
-9.999999999999999653057388335469287129029766072807834803722132479e-65
-0.003359382946410160745902574708741938902534229782759211957454681396
-9.999999999999999653057388335469287129029766072807834803722132479e-65
-9.999999999999999653057388335469287129029766072807834803722132479e-65
+1.871929361865491613930328508585976265976569266058504581451416015625e-06
+0.00023603394037677709348978238867999834837974049150943756103515625
+0.0037970706619493044438329309997470772941596806049346923828125
+0.050177628895479005344260059473526780493557453155517578125
+0.1098932324552358286151587662970996461808681488037109375
+0.007641542007296760717627304160259882337413728237152099609375
+0.0160148049969302298800588602034622454084455966949462890625
+0.19159429003178651829131240447168238461017608642578125
+0.177914330992157199151648683255189098417758941650390625
+0.38100416015747595732676700208685360848903656005859375
+7.502908426371362382163959607594971856769916485063731670379638671875e-07
+9.2212326230535645210824313711128752402146346867084503173828125e-05
+0.0003975693610413832197909844712313542913761921226978302001953125
+0.00550446768066507365058104284116780036129057407379150390625
+0.0293037459447654143052108821621004608459770679473876953125
+0.1998636693728008395520845397186349146068096160888671875
+0.006242539982388668751145388569057104177772998809814453125
+0.41671961732306883607890313214738853275775909423828125
+0.07706609115986708113954506416121148504316806793212890625
+0.26480933655832938899976625179988332092761993408203125
+1.0000000000000000540140885956810330905283414542659480243086298659334589074477275146251465374054496583672623295800994620910810169e-128
+1.0000000000000000540140885956810330905283414542659480243086298659334589074477275146251465374054496583672623295800994620910810169e-128
+0.0171191434004156380066508091886134934611618518829345703125
+0.24217601602121507919918030893313698470592498779296875
+0.59818906159696150215410170858376659452915191650390625
+0.2240890811920234482190750213703722693026065826416015625
+1.0000000000000000540140885956810330905283414542659480243086298659334589074477275146251465374054496583672623295800994620910810169e-128
+0.0033593829464101600430769867244862325605936348438262939453125
+1.0000000000000000540140885956810330905283414542659480243086298659334589074477275146251465374054496583672623295800994620910810169e-128
+1.0000000000000000540140885956810330905283414542659480243086298659334589074477275146251465374054496583672623295800994620910810169e-128
--- trunk/koffice/krita/plugins/viewplugins/painterlymixer/ipg/compile #697376:697377
@@ -1 +1 @@
-g++ -o ipg ipg.cpp -lm -lgmp -I/home/emanuele/devel/inst/kde/include
+g++ -o ipg ipg.cpp -lm -lgmp -lgmpxx -I/home/emanuele/devel/inst/kde/include
--- trunk/koffice/krita/plugins/viewplugins/painterlymixer/ipg/ipg.cpp #697376:697377
@@ -35,6 +35,7 @@
#include <gmm/gmm.h>
#include <gmp.h>
+#include <gmpxx.h>
#include "matching_curves.h"
@@ -44,73 +45,78 @@
#define PREC 1024
#define NUM 10
-class {
- long double mH[471];
- long double mW[NUM];
- long double mX[NUM];
- long double m_matrix[3][NUM];
+class c {
+ mpf_class mH[471];
+ mpf_class mW[NUM];
+ mpf_class mX[NUM];
+ mpf_class m_matrix[3][NUM];
public:
- long double &H(int i)
+ mpf_class &H(int i)
{
if (i >= 0 && i <= 470)
return mH[i];
}
- long double &W(int i)
+ mpf_class &W(int i)
{
if (i >= 0 && i < NUM)
return mW[i];
}
- long double &X(int i)
+ mpf_class &X(int i)
{
if (i >= 0 && i < NUM)
return mX[i];
}
- long double &matrix(int i, int j)
+ mpf_class &matrix(int i, int j)
{
if (i >= 0 && i < 3 && j >= 0 && j < NUM)
return m_matrix[i][j];
}
- const long double x(int i)
+ const mpf_class x(int i)
{
- if (i >= 0 && i <= 470)
- return cmf::x[i];
+ if (i >= 0 && i <= 470) {
+ mpf_class X(cmf::x[i]);
+ return X;
+ }
}
- const long double y(int i)
+ const mpf_class y(int i)
{
- if (i >= 0 && i <= 470)
- return cmf::y[i];
+ if (i >= 0 && i <= 470) {
+ mpf_class Y(cmf::y[i]);
+ return Y;
+ }
}
- const long double z(int i)
+ const mpf_class z(int i)
{
- if (i >= 0 && i <= 470)
- return cmf::z[i];
+ if (i >= 0 && i <= 470) {
+ mpf_class Z(cmf::z[i]);
+ return Z;
+ }
}
- const long double F(int i)
+ const mpf_class F(int i)
{
if (i >= 0 && i <= 470)
- return (x(i)+y(i)+z(i));
+ return x(i)+y(i)+z(i);
}
- const long double W_func(int i)
+ const mpf_class W_func(int i)
{
if (i >= 0 && i <= 470)
- return H(i) * F(i);
+ return H(i)*F(i);
}
} C;
-void init_v(mpf_t v[]);
-void clear_v(mpf_t v[]);
-void get_v(mpf_t v[], long double d[]);
-void pvalue(mpf_t v[], long double x, mpf_t res);
+void get_v(mpf_class v[], double d[]);
+void pvalue(mpf_class v[], mpf_class X, mpf_class &res);
-void copy(mpf_t from[], mpf_t to[]);
-void add(mpf_t v1[], mpf_t v2[], mpf_t vr[]);
-void sub(mpf_t v1[], mpf_t v2[], mpf_t vr[]);
-void mul(mpf_t v[], mpf_t k, mpf_t vr[]);
-void div(mpf_t v[], mpf_t k, mpf_t vr[]);
-void scalar(mpf_t v1[], mpf_t v2[], mpf_t r);
+void copy(mpf_class from[], mpf_class to[]);
+void add(mpf_class v1[], mpf_class v2[], mpf_class vr[]);
+void sub(mpf_class v1[], mpf_class v2[], mpf_class vr[]);
+void mul(mpf_class v[], mpf_class k, mpf_class vr[]);
+void div(mpf_class v[], mpf_class k, mpf_class vr[]);
+mpf_class pow(mpf_class d, int p);
+void scalar(mpf_class v1[], mpf_class v2[], mpf_class &r);
-long double standardIntegral(long double p);
+mpf_class standardIntegral(int p);
bool loadIlluminant(char *filename);
bool computeRoots();
bool computeWeights();
@@ -153,118 +159,102 @@
#define STD_FOR for (int i = 0; i <= NUM; i++)
-void init_v(mpf_t v[])
+void get_v(mpf_class v[], double d[])
{
STD_FOR
- mpf_init(v[i]);
+ d[i] = v[i].get_d();
}
-void clear_v(mpf_t v[])
+void pvalue(mpf_class v[], mpf_class X, mpf_class &res)
{
+ res = 0.0;
STD_FOR
- mpf_clear(v[i]);
+ res += pow(X,i) * v[i];
}
-void get_v(mpf_t v[], long double d[])
+void copy(mpf_class from[], mpf_class to[])
{
STD_FOR
- d[i] = mpf_get_d(v[i]);
+ to[i] = from[i];
}
-void pvalue(mpf_t v[], long double x, mpf_t res)
+void add(mpf_class v1[], mpf_class v2[], mpf_class vr[])
{
- mpf_t tmp, X;
-
- mpf_init(tmp);
- mpf_init(X);
-
- mpf_set_d(res, 0.0);
- STD_FOR {
- mpf_set_d(X, x);
- mpf_pow_ui(X, X, i);
- mpf_mul(tmp, v[i], X);
- mpf_add(res, res, tmp);
- }
-
- mpf_clear(tmp);
- mpf_clear(X);
-}
-
-void copy(mpf_t from[], mpf_t to[])
-{
STD_FOR
- mpf_set(to[i], from[i]);
+ vr[i] = v1[i] + v2[i];
}
-void add(mpf_t v1[], mpf_t v2[], mpf_t vr[])
+void sub(mpf_class v1[], mpf_class v2[], mpf_class vr[])
{
STD_FOR
- mpf_add(vr[i], v1[i], v2[i]);
+ vr[i] = v1[i] - v2[i];
}
-void sub(mpf_t v1[], mpf_t v2[], mpf_t vr[])
+void mul(mpf_class v[], mpf_class k, mpf_class vr[])
{
STD_FOR
- mpf_sub(vr[i], v1[i], v2[i]);
+ vr[i] = v[i] * k;
}
-void mul(mpf_t v[], mpf_t k, mpf_t vr[])
+void div(mpf_class v[], mpf_class k, mpf_class vr[])
{
STD_FOR
- mpf_mul(vr[i], v[i], k);
+ vr[i] = v[i] / k;
}
-void div(mpf_t v[], mpf_t k, mpf_t vr[])
+mpf_class pow(mpf_class d, int p)
{
- STD_FOR
- mpf_div(vr[i], v[i], k);
+ mpf_class res;
+
+ if (p == 0) {
+ res = 1;
+ return res;
+ }
+
+ res = d;
+
+ for (int i = 1; i < abs(p); i++)
+ res *= d;
+
+ if (p < 0)
+ res = 1.0/res;
+
+ return res;
}
-void scalar(mpf_t v1[], mpf_t v2[], mpf_t r)
+void scalar(mpf_class v1[], mpf_class v2[], mpf_class &r)
{
- mpf_t W, norm;
+ mpf_class W, norm;
- mpf_init(W);
- mpf_init(norm);
+ norm = 2.0/470.0;
- mpf_set_d(norm, 2.0/470.0);
-
- mpf_set_d(r, 0.0);
+ r = 0.0;
for (int i = 0; i <= 470; i++) {
- mpf_t val1, val2, curr;
- long double x = 2.0*(long double)(i)/470.0 - 1.0;
-// long double x = 360+i;
+ mpf_class val1, val2, x;
- mpf_init(val1);
- mpf_init(val2);
- mpf_init(curr);
+ x = 2.0*(double)(i)/470.0 - 1.0;
+ W = C.W_func(i);
- mpf_set_d(W, C.W_func(i));
pvalue(v1, x, val1);
pvalue(v2, x, val2);
- mpf_mul(curr, val1, val2);
- mpf_mul(curr, curr, W);
- mpf_mul(curr, curr, norm);
- mpf_add(r, r, curr);
-
- mpf_clear(curr);
- mpf_clear(val2);
- mpf_clear(val1);
+ r += W*val1*val2*norm;
}
-
- mpf_clear(norm);
- mpf_clear(W);
}
-long double standardIntegral(long double p)
+mpf_class standardIntegral(int p)
{
- const long double norm = 2.0/470.0;
- long double I = 0;
+ mpf_class norm, I;
+ norm = 2.0/470.0;
+ I = 0.0;
+
for (int i = 0; i <= 470; i++) {
- long double x = 2.0*(long double)(i)/470.0 - 1.0;
- I += C.W_func(i)*pow(x, p)*norm;
+ mpf_class tmp, x;
+
+ x = 2.0*(double)(i)/470.0 - 1.0;
+
+ I += pow(x, p)*C.W_func(i)*norm;
}
return I;
@@ -275,76 +265,61 @@
ifstream f; // TODO Add error handling
f.open(filename);
- for (int i = 0; i <= 470; i++)
- f >> C.H(i);
+ for (int i = 0; i <= 470; i++) {
+ double H;
+ f >> H;
+ C.H(i) = H;
+ }
return true;
}
bool computeRoots()
{
- mpf_t v[NUM+1][NUM+1];
- mpf_t u[NUM+1][NUM+1];
- mpf_t e[NUM+1][NUM+1];
+ mpf_class v[NUM+1][NUM+1];
+ mpf_class u[NUM+1][NUM+1];
+ mpf_class e[NUM+1][NUM+1];
- STD_FOR {
- init_v(v[i]);
- init_v(u[i]);
- init_v(e[i]);
- }
-
for (int i = 0; i <= NUM; i++) {
- mpf_set_d(v[i][i], 1.0);
+ v[i][i] = 1.0;
copy(v[i], u[i]);
cout << "v = [";
for (int k = 0; k < NUM; k++)
- cout << mpf_get_d(v[i][k]) << ",";
- cout << mpf_get_d(v[i][NUM]) << "]" << endl;
+ cout << v[i][k].get_d() << ",";
+ cout << v[i][NUM].get_d() << "]" << endl;
for (int j = 0; j < i; j++) {
- mpf_t et[NUM+1];
- mpf_t st;
+ mpf_class et[NUM+1];
+ mpf_class st;
- init_v(et);
- mpf_init(st);
-
scalar(e[j], v[i], st);
mul(e[j], st, et);
sub(u[i], et, u[i]);
-
- mpf_clear(st);
- clear_v(et);
-
}
cout << "u = [";
for (int k = 0; k < NUM; k++)
- cout << mpf_get_d(u[i][k]) << ",";
- cout << mpf_get_d(u[i][NUM]) << "]" << endl;
+ cout << u[i][k].get_d() << ",";
+ cout << u[i][NUM].get_d() << "]" << endl;
- mpf_t n;
+ mpf_class n;
- mpf_init(n);
-
scalar(u[i], u[i], n);
- mpf_sqrt(n, n);
+ n = sqrt(n);
div(u[i], n, e[i]);
cout << "e = [";
for (int k = 0; k < NUM; k++)
- cout << mpf_get_d(e[i][k]) << ",";
- cout << mpf_get_d(e[i][NUM]) << "]" << endl;
+ cout << e[i][k].get_d() << ",";
+ cout << e[i][NUM].get_d() << "]" << endl;
cout << "----------------------------------" << endl;
-
- mpf_clear(n);
-
}
- long double d[NUM+1];
+ double d[NUM+1];
get_v(e[NUM], d);
@@ -365,21 +340,15 @@
cout << "IN [-1, 1]: ";
for (int i = 0; i < NUM; i++) {
- long double X;
+ double X;
in >> X;
C.X(i) = X;
- cout << C.X(i) << ", ";
+ cout << X << ", ";
}
cout << endl;
in.close();
- STD_FOR {
- clear_v(v[i]);
- clear_v(u[i]);
- clear_v(e[i]);
- }
-
system("rm /tmp/ipg.tmp");
system("rm /tmp/ipg_ans_raw.tmp");
system("rm /tmp/ipg_ans.tmp");
@@ -389,13 +358,13 @@
bool computeWeights()
{
- dense_matrix<long double> M(NUM, NUM);
- vector<long double> B(NUM), W(NUM);
+ dense_matrix<mpf_class> M(NUM, NUM);
+ vector<mpf_class> B(NUM), W(NUM);
clear(M);
for (int i = 0; i < NUM; i++) {
- B[i] = standardIntegral((long double)i);
+ B[i] = standardIntegral(i);
for (int j = 0; j < NUM; j++)
- M(i, j) = pow(C.X(j), (long double)i);
+ M(i,j) = pow(C.X(j), i);
}
lu_solve(M, W, B);
cout << B << endl << M << endl;
@@ -404,8 +373,9 @@
for (int i = 0; i < NUM; i++) {
C.W(i) = W[i];
// Integer-ize abscissas
- C.X(i) = (int)(470.0*C.X(i)/2.0 + 470.0/2.0);
- cout << C.X(i) << ", ";
+ int X = (int)(470.0*C.X(i).get_d()/2.0 + 470.0/2.0);
+ C.X(i) = X;
+ cout << X << ", ";
}
cout << endl;
@@ -415,10 +385,10 @@
bool computeMatrix()
{
cout.precision(128);
- long double Hy_integral = 0, k;
+ mpf_class Hy_integral(0.0), k;
for (int i = 0; i < NUM; i++) {
- int X = (int)C.X(i);
+ int X = C.X(i).get_ui();
Hy_integral += C.W(i)*C.y(X)/C.F(X);
}
@@ -428,18 +398,18 @@
cout << "K: " << k << endl;
for (int j = 0; j < NUM; j++) {
- int X = (int)C.X(j);
+ int X = C.X(j).get_ui();
C.matrix(0,j) = k*(470.0/2.0)*C.W(j)*C.x(X)/C.F(X);
}
for (int j = 0; j < NUM; j++) {
- int X = (int)C.X(j);
+ int X = C.X(j).get_ui();
C.matrix(1,j) = k*(470.0/2.0)*C.W(j)*C.y(X)/C.F(X);
}
for (int j = 0; j < NUM; j++) {
- int X = (int)C.X(j);
+ int X = C.X(j).get_ui();
C.matrix(2,j) = k*(470.0/2.0)*C.W(j)*C.z(X)/C.F(X);
}
-
+/*
cout.precision(128);
cout << "double g_matrix[3][10] = {" << endl;
for (int i = 0; i < 3; i++) {
@@ -456,7 +426,7 @@
cout << endl;
}
cout << "};" << endl;
-
+*/
return true;
}
@@ -466,13 +436,13 @@
f.open(filename);
- f.precision(64);
+ f.precision(128);
for (int i = 0; i < 3; i++)
- for (int j = 0; j < 10; j++)
- if (C.matrix(i,j))
- f << C.matrix(i,j) << endl;
+ for (int j = 0; j < NUM; j++)
+ if (C.matrix(i,j).get_d())
+ f << C.matrix(i,j).get_d() << endl;
else
- f << 1e-64 << endl;
+ f << 1e-128 << endl;
return true;
}
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic