Programowanie w systemie UNIX/qd

Nazwy

• pakiety
  • double-double
  • quad-double package
• biblioteka : qd-library
• precyzja
  • double- double precision
  • quaddruple-precision

API

• c
• c++
• fortran

 "Most of the routines runs noticeably faster if implemented in C" then in C++

wrappers

• python

typy liczb (precyzja)

• quad-double = qd_real = four times the precision of IEEE double = 212 mantissa bits = approximately 64 decimal digits
• double-double = dd_real = twice the precision of IEEE double = 106 mantissa bits = approximately 32 decimal digits
• double

// from qd_inline.h
inline double to_double(const qd_real &a) {
  return a[0];
}

instalacja

pliki

• libqd.a
• qd.h

zależności:

• autoconf
• m4
• git

git

git clone git@bitbucket.org:carter/libqd-reference.git
 autoreconf --install
 autoconf
 ./configure
 make 
 sudo make install
 

Biblioteka została zainstalowana do :

  /usr/local/lib

Sprawdzamy :

 ldconfig -p | grep libqd

przykładowy wynik :

 libqd.so.0 (libc6,x86-64) => /usr/local/lib/libqd.so.0
 libqd.so (libc6,x86-64) => /usr/local/lib/libqd.so

pakiet

• libqd
• libqd-dev

 sudo apt-get update
 sudo apt-get install libqd-dev

funkcje

znak liczby

 #include <qd/c_qd.h> // qd library :  quad double number with precision = 64 decimal digits 

 double qd[4]; // qd number   
 // compute qd 
 if( qd[0]>0){}; // check sign of double(qd)

dodawanie

 c_qd_add(a, b, c); // Put (a+b) into c ( c = a+b )

write

Jak wy świetlić liczbę qd  ?

 #include <qd/c_qd.h> // qd library :  quad double number with precision = 64 decimal digits 

 double qd[4]; // qd number   
 // compute qd 
 // print qd
 printf("qd = "); c_qd_write(qd); // 64 decimal digits
 printf("qd = %.16f\n", qd[0]);  // 16 decimal digits

Jak zapisać liczbę qd do pliku ?

#include <qd/c_qd.h> // qd library :  quad double number with precision = 64 decimal digits 


 double qd[4]; // qd number  
 char   qds[1000]; // string for qd number

 c_qd_swrite (qd, 64, qds, 1000); // copy d to string ds
 fprintf(file, " %s \n", qds);  // save ds to the file

Użycie

Obraz utworzony z pomocą programu używajaćego libqd

compileExample.cpp

// simple exaple of QD usage to illustrate linking process
// Alex Kaiser, LBNL, 6/3/2010



#include <iostream>
#include <qd/qd_real.h>
#include <qd/fpu.h>

using namespace std; 

int main() {
	
	// ensure that 80-bit arithmetic is not in place
	// this call forces 64-bit arithmetic
	unsigned int old_cw;
	fpu_fix_start(&old_cw);

	cout.precision(60); 
	
	// simple read example
	/*
	qd_real readTest ; 
	cin >> readTest ; 	
	cout << "readTest = " << readTest << endl ;
	 */ 
	
	// simple demo
	qd_real x = "1.0" ;  
	x /= 3.0 ;
	qd_real y ; 
	y = pow( qd_real(2.0) , 3 ) ; 
	cout << "y = " << y << endl;
	cout << "x = " << x << endl; 
	
	
	qd_real a ; 
	qd_real b = qd_real("0.1");
	
	a = sqrt(b);
	cout << " sqrt(0.1) = " << a << endl;
	cout << " sqrt(0.1) * sqrt(0.1) = " << a * a << endl; 
	
	fpu_fix_end(&old_cw); 
	return 0;
}

 cd tests
 g++ compileExample.cpp  -lqd
 ./a.out
 y = 8.000000000000000000000000000000000000000000000000000000000000e+00
 x = 3.333333333333333333333333333333333333333333333333333333333333e-01
 sqrt(0.1) = 3.162277660168379331998893544432718533719555139325216826857505e-01
 sqrt(0.1) * sqrt(0.1) = 1.000000000000000000000000000000000000000000000000000000000000e-01

pi

/*
 gcc p.c -lqd -Wall
 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib/
 ./a.out 


*/

#include <stdio.h>
#include <math.h>
#include <qd/c_qd.h>



/* Test 1.  Salamin-Brent quadratically convergent formula for pi. */
int test_1() {

  double a[4], b[4], s[4], p[4], t[4], t2[4];
  double a_new[4], b_new[4], p_old[4];
  double m, err;
  int r, i;
  const int max_iter = 20;

  puts("Test 1.  (Salamin-Brent quadratically convergent formula for pi)");

  c_qd_copy_d(1.0, a);  /* a = 1.0 */
  c_qd_copy_d(0.5, t);  /* t = 0.5 */
  c_qd_sqrt(t, b);      /* b = sqrt(t) */
  c_qd_copy_d(0.5, s);  /* s = 0.5 */
  m = 1.0;

  c_qd_sqr(a, p);
  c_qd_selfmul_d(2.0, p);
  c_qd_selfdiv(s, p);

  printf("  iteration 0: ");
  c_qd_write(p);
  for (i = 1; i <= max_iter; i++) {
    m *= 2.0;

    /* a_new = 0.5 * (a + b) */
    c_qd_add(a, b, a_new);
    c_qd_selfmul_d(0.5, a_new);

    c_qd_mul(a, b, b_new); /* b_new = a * b */

    /* Compute s = s - m * (a_new^2 - b) */
    c_qd_sqr(a_new, t);       /* t = a_new ^ 2 */
    c_qd_selfsub(b_new, t);   /* t -= b_new */
    c_qd_selfmul_d(m, t);     /* t *= m */
    c_qd_selfsub(t, s);       /* s -= t */

    c_qd_copy(a_new, a);
    c_qd_sqrt(b_new, b);
    c_qd_copy(p, p_old);

    /* Compute  p = 2.0 * a^2 / s */
    c_qd_sqr(a, p);
    c_qd_selfmul_d(2.0, p);
    c_qd_selfdiv(s, p);

    /* Test for convergence by looking at |p - p_old|. */
    c_qd_sub(p, p_old, t);
    c_qd_abs(t, t2);
    c_qd_comp_qd_d(t2, 1e-60, &r);
    if (r < 0) break;

    printf("  iteration %1d: ", i);
    c_qd_write(p);
  }

  c_qd_pi(p);   /* p = pi */
  printf("          _pi: ");
  c_qd_write(p);
  printf("        error: %.5e = %g eps\n", t2[0], t2[0] / ldexp(1.0, -209));

  return 0;
}





int main(void) {
  fpu_fix_start(NULL);
  return test_1();
}

Wynik:

./a.out
Test 1.  (Salamin-Brent quadratically convergent formula for pi)
  iteration 0: 4.00000000000000000000000000000000000000000000000000000000000000e+00
  iteration 1: 3.18767264271210862720192997052536923265105357185936922648763399e+00
  iteration 2: 3.14168029329765329391807042456000938279571943881540283264418946e+00
  iteration 3: 3.14159265389544649600291475881804348610887923726131158965110136e+00
  iteration 4: 3.14159265358979323846636060270663132175770241134242935648684602e+00
  iteration 5: 3.14159265358979323846264338327950288419716994916472660583469612e+00
  iteration 6: 3.14159265358979323846264338327950288419716939937510582097494458e+00
          _pi: 3.14159265358979323846264338327950288419716939937510582097494459e+00
        error: 8.20417e-63 = 6.75 eps

Złota proporcja

/*

gcc c.c -lqd -Wall -DUSE_QD_REAL
 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib/
 ./a.out 



1.61803398874989484820458683436563811772030917980576286213544862e+00


*/

#include <stdio.h>
#include <math.h>
#include <qd/c_qd.h> // 64 decimal digits 

// phi = (1+sqrt(5))/2 = golden ratio
// https://en.wikipedia.org/wiki/Golden_ratio
int GivePhi(double p[4]){
 

  double a[4];
  
  


  c_qd_copy_d(5.0, p); // a = 5.0
  c_qd_copy_d(1.0, a); // b = 1.0

  c_qd_sqrt(p, p); // a = sqrt(a) =sqrt(5)
  c_qd_add(p, a, p); // a= a+b = 1+sqrt(5)
  c_qd_selfmul_d(0.5, p); // a = a*0.5 = a/2 = phi
  
  return 0;

}





int main(void) {

  double phi[4];
  
  fpu_fix_start(NULL); // turns on the round-to-double bit in the FPU control word on Intel x86 Processors. 

  GivePhi(phi); 
  //
  printf("phi = ");
  c_qd_write(phi); // print 
  printf("\n");



  //
  fpu_fix_end(NULL);

  return 0; 
}

ułamki łańcuchowe

/*

gcc c.c -lqd -Wall -DUSE_QD_REAL
 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib/
 ./a.out 



1.61803398874989484820458683436563811772030917980576286213544862e+00
./a.out
t(0)  = 2.95685999407889202537083517598191479880016272621355940112392033e-01
t(1)  = 2.87561745861029587995763349374215634699576366286473943622953355e-01
t(2)  = 2.85916253540726005296290351777281930430489555597203400050155169e-01
t(3)  = 2.85734672540588170268886130620030074831025799037914834916050614e-01
t(4)  = 2.85716326317041655001121670485871240052920659174284332597727197e-01
t(5)  = 2.85714489793782460278353122604001584582831972498394265344044955e-01
t(6)  = 2.85714306122427620319960416932943500595980573209849938259909447e-01
t(7)  = 2.85714287755101827223406574888790339883583192331314008460721020e-01
t(8)  = 2.85714285918367344802846109454092778340593443209054635310671634e-01
t(9)  = 2.85714285734693877529661113954572642424219379874139361821960174e-01
t(10) = 2.85714285716326530612031305016895554055165716643985018539739541e-01
t(11) = 2.85714285714489795918365211009352427817161964062263710683311655e-01
2/7   = 2.85714285714285714285714285714285714285714285714285714285714286e-01


*/

#include <stdio.h>
#include <math.h>
#include <qd/c_qd.h> // 64 decimal digits 



// phi = (1+sqrt(5))/2 = golden ratio
// https://en.wikipedia.org/wiki/Golden_ratio
int GivePhi(double p[4]){
 

  double a[4];
  
  


  c_qd_copy_d(5.0, p); // p = 5.0
  c_qd_copy_d(1.0, a); // a = 1.0

  c_qd_sqrt(p, p); // p = sqrt(p) =sqrt(5)
  c_qd_add(p, a, p); // p= a+p = 1+sqrt(5)
  c_qd_selfmul_d(0.5, p); // p = p*0.5 = p/2 = phi
  
  return 0;

}



/*
 Real tt(int n)
{
  // t(n) = [0;3, 2 , 10^n, golden_ratio]
  Real phi = (sqrt(Real(5)) + 1) / 2;
  Real tenn = pow(Real(10), n);
  return 0 +1/( 3 +1/( 2 +1/( tenn +1/( phi ))));
}


*/



void GiveT(int n,  double t[4]){

double p[4], a[4], b[4], one[4], two[4], three[4];

double phi[4];

  

 c_qd_copy_d(1.0, one); // one = 1.0
 //c_qd_write(one); // print 
 c_qd_copy_d(2.0, two); // one = 1.0
 //c_qd_write(two); // print 
 c_qd_copy_d(3.0, three); // one = 1.0
 //c_qd_write(three); // print 
 //printf("\n");

 //c_qd_selfdiv(three,two); //two=two/three
 //c_qd_write(two); // print 

 //printf("10 \n");
 c_qd_copy_d(10.0, a); // a = 10.0
 //c_qd_write(a); // print  

 //printf("10^n \n");
 c_qd_npwr(a, n, p); // p = a^n
 //c_qd_write(p); // print 


 GivePhi(phi); 
 //printf("phi = ");
 //c_qd_write(phi); // print 
 //printf("\n");


 //printf("1/phi \n");
 c_qd_div(one,phi,b); // b= one/phi = 1/phi
 //c_qd_write(b); // print 

 //printf("10^n + 1/phi \n");
 c_qd_add(p,b,t); // t= p+b = p+1/phi = (a^n) +(1/phi)
 //c_qd_write(t); // print 

 //printf("1/(10^n + 1/phi) \n");
 c_qd_copy(one,b);
 c_qd_selfdiv(t,b); //b= 1/t
 //c_qd_write(b); // print 

 //printf("2+(1/(10^n + 1/phi)) \n");
 c_qd_add(two,b,t);
 //c_qd_write(t); // print 

  //printf("1/(2+(1/(10^n + 1/phi))) \n");
 c_qd_copy(one,b);
 c_qd_selfdiv(t,b); // 

 //printf("3+1/(2+(1/(10^n + 1/phi))) \n");
 c_qd_add(three,b, t); //  
// c_qd_write(t); // print 
 

 //printf("1/(3+1/(2+(1/(10^n + 1/phi)))) \n");
 c_qd_copy(one,b);
 c_qd_selfdiv(t,b); //t= 1/b
 //c_qd_write(b); // print 

 c_qd_copy(b,t);
  
 

}


void Rat2float(double numerator, double denominator, double result[4]){

 double n[4];
 double d[4];

 c_qd_copy_d(numerator, n); // one = 1.0
 c_qd_copy_d(denominator, d); // one = 1.0

 c_qd_div(n, d, result); // b= one/phi = 1/phi
 c_qd_write(result); // print 


}




int main(void) {

  
  double t[4];
  int n;
  
  fpu_fix_start(NULL); // turns on the round-to-double bit in the FPU control word on Intel x86 Processors. 

  
  for (n=0; n<12; n++){
  
  GiveT(n,t);
  printf("t(%d) = ", n);
  c_qd_write(t); // print 
  printf("\n");
}
  Rat2float( (double) 2, (double) 7, t);  
  
  //
  fpu_fix_end(NULL);

  return 0; 
}

Wynik :

t(0)  = 2.95685999407889202537083517598191479880016272621355940112392033e-01
t(1)  = 2.87561745861029587995763349374215634699576366286473943622953355e-01
t(2)  = 2.85916253540726005296290351777281930430489555597203400050155169e-01
t(3)  = 2.85734672540588170268886130620030074831025799037914834916050614e-01
t(4)  = 2.85716326317041655001121670485871240052920659174284332597727197e-01
t(5)  = 2.85714489793782460278353122604001584582831972498394265344044955e-01
t(6)  = 2.85714306122427620319960416932943500595980573209849938259909447e-01
t(7)  = 2.85714287755101827223406574888790339883583192331314008460721020e-01
t(8)  = 2.85714285918367344802846109454092778340593443209054635310671634e-01
t(9)  = 2.85714285734693877529661113954572642424219379874139361821960174e-01
t(10) = 2.85714285716326530612031305016895554055165716643985018539739541e-01
t(11) = 2.85714285714489795918365211009352427817161964062263710683311655e-01
2/7   = 2.85714285714285714285714285714285714285714285714285714285714286e-01

inicjalizacja

W c++

• strings must be used to assign to a quad-double (or double-double) numbers
• otherwise the double precision approximation is assigned.

  qd_real a = "3.141592653589793238462643383279502884197169399375105820";
  dd_real b = "2.249775724709369995957";
  double c = 0.1

stałe

 qd real:: pi
 qd real:: pi2
 qd real:: pi4
 qd real:: e
 qd real:: log2

  // qd_real.h
  static const qd_real _2pi;
  static const qd_real _pi;
  static const qd_real _3pi4;
  static const qd_real _pi2;
  static const qd_real _pi4;
  static const qd_real _e;
  static const qd_real _log2;
  static const qd_real _log10;
  static const qd_real _nan;
  static const qd_real _inf;

Pomoc

• David Bailey's double-double package
• Yozo Hida library
• https://www.researchgate.net/publication/228570156_Library_for_Double-Double_and_Quad-Double_Arithmetic
• https://www.codeproject.com/articles/884606/the-double-double-type
• http://web.mit.edu/tabbott/Public/quaddouble-debian/qd-2.3.4-old/docs/qd.pdf
• kod żródłowy

Zobacz również

• libquadmath - GNU C extension
• gpuprec = QD na GPU.
  • https://code.google.com/archive/p/gpuprec/
  • http://homepages.math.uic.edu/~jan/mcs572/mcs572notes/lec37.html
  • http://homepages.math.uic.edu/~jan/mcs572/quad_double_cuda.pdf
  • https://github.com/lumianph/gpuprec