Programowanie w systemie UNIX/MPFR
Wymagania
edytuj- GMP
- dla tworzenia dokumentacji
- texi2dvi (texinfo)
- texlive
Instalacja
edytujPo zainstalowaniu GMP instalujemy MPFR za pomocą:
- graficznego menadżera pakietów Synaptic w Ubuntu i Debianie
- apt-get (
sudo apt-get install libmpfr-dev libmpfr-doc libmpfr4 libmpfr4-dbg) - ze źródeł
Opis instalacji ze źródeł jest w pliku install:
- skopiuj aktualny plik archiwum (np.:
mpfr-3.1.2.tar.xz) - rozpakuj plik
- zastosuj łaty
- zbadaj system (
./configure) - zbuduj bibliotekę (
make) - zrób testy
- zainstaluj
curl http://www.mpfr.org/mpfr-3.1.2/allpatches | patch -N -Z -p1 ./configure make make check sudo make install make pdf
Położenie plików / katalogi
edytujPołożenie bibliotek
edytujsudo updatedb locate libmpfr
przykładowy wynik:
/home/a/Pobrane/mpfr-3.1.5/src/libmpfr.la /home/a/Pobrane/mpfr-3.1.5/src/.libs/libmpfr.a /home/a/Pobrane/mpfr-3.1.5/src/.libs/libmpfr.la /home/a/Pobrane/mpfr-3.1.5/src/.libs/libmpfr.lai /home/a/Pobrane/mpfr-3.1.5/src/.libs/libmpfr.so /home/a/Pobrane/mpfr-3.1.5/src/.libs/libmpfr.so.4 /home/a/Pobrane/mpfr-3.1.5/src/.libs/libmpfr.so.4.1.5 /usr/lib/x86_64-linux-gnu/libmpfr.so.4 /usr/lib/x86_64-linux-gnu/libmpfr.so.4.1.4 /usr/local/lib/libmpfr.a /usr/local/lib/libmpfr.la /usr/local/lib/libmpfr.so /usr/local/lib/libmpfr.so.4 /usr/local/lib/libmpfr.so.4.1.5 /usr/share/doc/libmpfr4 /usr/share/doc/libmpfr4/AUTHORS /usr/share/doc/libmpfr4/BUGS /usr/share/doc/libmpfr4/NEWS.gz /usr/share/doc/libmpfr4/README /usr/share/doc/libmpfr4/README.Debian /usr/share/doc/libmpfr4/TODO.gz /usr/share/doc/libmpfr4/changelog.Debian.gz /usr/share/doc/libmpfr4/copyright /var/lib/dpkg/info/libmpfr4:amd64.list /var/lib/dpkg/info/libmpfr4:amd64.md5sums /var/lib/dpkg/info/libmpfr4:amd64.shlibs /var/lib/dpkg/info/libmpfr4:amd64.symbols
Pliki biblioteki są standardowo instalowane w:
- /usr/local/lib/
- libmpfr.a ( biblioteka statyczna)
- libmpfr.so ( biblioteka dynamiczna )
- /usr/lib/x86_64-linux-gnu ( dla architektury 64 bitowej, zobacz MultiArchSpec)
Wersje pakietów i bibliotek
- libmpfr6 zawiera wersję MPFR 4.0.1-1
- libmpfr4 zzwiera wersję MPFR 3.1.2-1
Sprawdzenie :
ldconfig -p | grep mpfr
Przykładowy wynik:
libmpfr.so.6 (libc6,x86-64) => /usr/lib/x86_64-linux-gnu/libmpfr.so.6 libmpfr.so.4 (libc6,x86-64) => /usr/local/lib/libmpfr.so.4 libmpfr.so (libc6,x86-64) => /usr/local/lib/libmpfr.so libmpfr.so (libc6,x86-64) => /usr/lib/x86_64-linux-gnu/libmpfr.so
Pliki nagłówkowe
edytujPlik nagłówkowy mpfr.h powinien być w /usr/include
/*
gcc d.c -lmpfr -lgmp
https://tspiteri.gitlab.io/gmp-mpfr-sys/mpfr/Installing-MPFR.html
gcc d.c -lmpfr -lgmp -Wall
./a.out
MPFR library: 4.0.1
MPFR header: 3.1.5 (based on 3.1.5)
*/
#include <stdio.h>
#include <mpfr.h>
int main (void)
{
printf ("MPFR library: %-12s\nMPFR header: %s (based on %d.%d.%d)\n",
mpfr_get_version (), MPFR_VERSION_STRING, MPFR_VERSION_MAJOR,
MPFR_VERSION_MINOR, MPFR_VERSION_PATCHLEVEL);
gmp_printf (" GMP-%s \n ", gmp_version );
return 0;
}
Problemy
edytuj- sudo apt-get remove libmpfr4
- remove old so files
wartości specjalne
edytujwartości specjalne (ang. special values)[1]:
- signed zeros = +0, -0
- infinities =
- not-a-number = NaN
typy
edytuj- mpf_t = typ z GMP dla liczb zmiennoprzecinkowych
- mpfr_t = float ("A floating-point number, or float for short, is an arbitrary precision significand (also called mantissa) with a limited precision exponent" [2])
- mpfr_prec_t =
- mpfr_rnd_t =
- ui = unsigned long int
- si = signed long int
- d = double
- ld = long double
- typy GMP dla liczb wielokrotnej precyzji:
- z oznacza mpz_t dla liczb całkowitych (ang. a multiple precision integer)
- q oznacza mpq_t dla liczby wymiernych (c nie ma takiego typu)
- f oznacza mpf_t dla liczb zmiennoprzecinkowych (ang. float = an arbitrary precision mantissa with a limited precision exponent)
typedef __mpfr_struct mpfr_t[1];
precision
edytuj- liczba bitów przeznaczonych do przedstawiania mantyssy ( ang. sigificand or mantissa)
- MPFR_PREC_MIN <= integer <= MPFR_PREC_MAX
funkcje
edytujmpfr_zero_p
edytujZnajdujemy definicję:
grep -nR "mpfr_zero_p"
w wyniku znajdujemy:
src/mpfr.h:753:#define mpfr_zero_p(_x) ((_x)->_mpfr_exp == __MPFR_EXP_ZERO)
czyli definicja:
#define mpfr_zero_p(_x) ((_x)->_mpfr_exp == __MPFR_EXP_ZERO) #define __MPFR_EXP_MAX ((mpfr_exp_t) (((mpfr_uexp_t) -1) >> 1)) #define __MPFR_EXP_ZERO (0 - __MPFR_EXP_MAX)
schemat programu
edytuj- deklaracja
- inicjalizacja
- nadanie wartości
- obliczenia
- czyszczenie pamięci
Pierwszy program
edytujonline
edytujMożna wypróbować MPFR online[3]
wersja biblioteki
edytujmpfr_printf(" MPFR-%s \n GMP-%s \n", mpfr_version, gmp_version );
granica ciągu
edytujProgram oblicza dolną granicę ciągu: 1+1/1!+1/2!+...+1/100!
używając 200-bit precyzji [4]
#include <stdio.h>
#include <gmp.h>
#include <mpfr.h>
int main (void)
{
unsigned int i;
mpfr_t s, t, u;
mpfr_init2 (t, 200);
mpfr_set_d (t, 1.0, MPFR_RNDD);
mpfr_init2 (s, 200);
mpfr_set_d (s, 1.0, MPFR_RNDD);
mpfr_init2 (u, 200);
for (i = 1; i <= 100; i++)
{
mpfr_mul_ui (t, t, i, MPFR_RNDU);
mpfr_set_d (u, 1.0, MPFR_RNDD);
mpfr_div (u, u, t, MPFR_RNDD);
mpfr_add (s, s, u, MPFR_RNDD);
}
printf ("Sum is ");
mpfr_out_str (stdout, 10, 0, s, MPFR_RNDD);
putchar ('\n');
mpfr_clear (s);
mpfr_clear (t);
mpfr_clear (u);
return 0;
}
Zapisujemy jako s.c i kompilujemy:
gcc s.c -lmpfr -lgmp
Uruchamiamy:
./a.out
i otrzymujemy wynik:
Sum is 2.7182818284590452353602874713526624977572470936999595749669131
Więcej
edytuj- użycie biblioteki MPFR i metody Newtona obliczeń i rysowania zbioru Mandelbrota [5]
znajdowanie punktu centralnego składowej zbioru Mandelbrota
edytuj/*
code from unnamed c program ( book )
http://code.mathr.co.uk/book
see mandelbrot_nucleus.c
by Claude Heiland-Allen
COMPILE :
gcc -std=c99 -Wall -Wextra -pedantic -O3 -ggdb m.c -lm -lmpfr
usage:
./a.out bits cx cy period maxiters
output : space separated complex nucleus on stdout
example
./a.out 53 -1.75 0 3 100
./a.out 53 -1.75 0 30 100
./a.out 53 0.471 -0.3541 14 100
./a.out 53 -0.12 -0.74 3 100
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h> // arbitrary precision
#include <mpfr.h>
extern int mandelbrot_nucleus(mpfr_t cx, mpfr_t cy, const mpfr_t c0x, const mpfr_t c0y, int period, int maxiters) {
int retval = 0;
mpfr_t zx, zy, dx, dy, s, t, u, v;
mpfr_inits2(mpfr_get_prec(c0x), zx, zy, dx, dy, s, t, u, v, (mpfr_ptr) 0);
mpfr_set(cx, c0x, GMP_RNDN);
mpfr_set(cy, c0y, GMP_RNDN);
for (int i = 0; i < maxiters; ++i) {
// z = 0
mpfr_set_ui(zx, 0, GMP_RNDN);
mpfr_set_ui(zy, 0, GMP_RNDN);
// d = 0
mpfr_set_ui(dx, 0, GMP_RNDN);
mpfr_set_ui(dy, 0, GMP_RNDN);
for (int p = 0; p < period; ++p) {
// d = 2 * z * d + 1;
mpfr_mul(u, zx, dx, GMP_RNDN);
mpfr_mul(v, zy, dy, GMP_RNDN);
mpfr_sub(u, u, v, GMP_RNDN);
mpfr_mul_2ui(u, u, 1, GMP_RNDN);
mpfr_mul(dx, dx, zy, GMP_RNDN);
mpfr_mul(dy, dy, zx, GMP_RNDN);
mpfr_add(dy, dx, dy, GMP_RNDN);
mpfr_mul_2ui(dy, dy, 1, GMP_RNDN);
mpfr_add_ui(dx, u, 1, GMP_RNDN);
// z = z^2 + c;
mpfr_sqr(u, zx, GMP_RNDN);
mpfr_sqr(v, zy, GMP_RNDN);
mpfr_mul(zy, zx, zy, GMP_RNDN);
mpfr_sub(zx, u, v, GMP_RNDN);
mpfr_mul_2ui(zy, zy, 1, GMP_RNDN);
mpfr_add(zx, zx, cx, GMP_RNDN);
mpfr_add(zy, zy, cy, GMP_RNDN);
}
// check d == 0
if (mpfr_zero_p(dx) && mpfr_zero_p(dy)) {
retval = 1;
goto done;
}
// st = c - z / d
mpfr_sqr(u, dx, GMP_RNDN);
mpfr_sqr(v, dy, GMP_RNDN);
mpfr_add(u, u, v, GMP_RNDN);
mpfr_mul(s, zx, dx, GMP_RNDN);
mpfr_mul(t, zy, dy, GMP_RNDN);
mpfr_add(v, s, t, GMP_RNDN);
mpfr_div(v, v, u, GMP_RNDN);
mpfr_mul(s, zy, dx, GMP_RNDN);
mpfr_mul(t, zx, dy, GMP_RNDN);
mpfr_sub(zy, s, t, GMP_RNDN);
mpfr_div(zy, zy, u, GMP_RNDN);
mpfr_sub(s, cx, v, GMP_RNDN);
mpfr_sub(t, cy, zy, GMP_RNDN);
// uv = st - c
mpfr_sub(u, s, cx, GMP_RNDN);
mpfr_sub(v, t, cy, GMP_RNDN);
// c = st
mpfr_set(cx, s, GMP_RNDN);
mpfr_set(cy, t, GMP_RNDN);
// check uv = 0
if (mpfr_zero_p(u) && mpfr_zero_p(v)) {
retval = 2;
goto done;
}
} // for (int i = 0; i < maxiters; ++i)
done: mpfr_clears(zx, zy, dx, dy, s, t, u, v, (mpfr_ptr) 0);
return retval;
}
void DescribeStop(int stop)
{
switch( stop )
{
case 0:
printf(" method stopped because i = maxiters\n");
break;
case 1:
printf(" method stopped because derivative == 0\n");
break;
//...
case 2:
printf(" method stopped because uv = 0\n");
break;
}
}
void usage(const char *progname) {
fprintf(stderr,
"program finds one center ( nucleus) of hyperbolic component of Mandelbrot set using Newton method"
"usage: %s bits cx cy period maxiters\n"
"outputs space separated complex nucleus on stdout\n"
"example %s 53 -1.75 0 3 100\n",
progname, progname);
}
int main(int argc, char **argv) {
// check the input
if (argc != 6) { usage(argv[0]); return 1; }
// read the values
int bits = atoi(argv[1]);
mpfr_t cx, cy, c0x, c0y;
mpfr_inits2(bits, cx, cy, c0x, c0y, (mpfr_ptr) 0);
mpfr_set_str(c0x, argv[2], 10, GMP_RNDN);
mpfr_set_str(c0y, argv[3], 10, GMP_RNDN);
int period = atoi(argv[4]);
int maxiters = atoi(argv[5]);
int stop;
//
stop = mandelbrot_nucleus(cx, cy, c0x, c0y, period, maxiters);
//
printf(" nucleus ( center) of component with period = %s near c = %s ; %s is : \n ", argv[4], argv[2], argv[3]);
mpfr_out_str(0, 10, 0, cx, GMP_RNDN);
putchar(' ');
mpfr_out_str(0, 10, 0, cy, GMP_RNDN);
putchar('\n');
//
DescribeStop(stop) ;
// clear memeory
mpfr_clears(cx, cy, c0x, c0y, (mpfr_ptr) 0);
//
return 0;
}
Znajdowanie kątów zewnętrznych promieni które lądują na korzeniach głównej składowej zbioru Mandelbrota
edytuj/*
------- Git -----------------
cd existing_folder
git init
git remote add origin git@gitlab.com:adammajewski/wake_gmp.git
git add .
git commit -m ""
git push -u origin master
-------------------------------
?? http://stackoverflow.com/questions/2380415/how-to-cut-a-mpz-t-into-two-parts-using-gmp-lib-on-c
to compile from console:
gcc w.c -lgmp -lmpfr -Wall
to run from console :
./a.out
tested on Ubuntu 14.04 LTS
uiIADenominator = 89
Using MPFR-3.1.2-p3 with GMP-5.1.3 with precision = 200 bits
internal angle = 34/89
first external angle :
period = denominator of internal angle = 89
external angle as a decimal fraction = 179622968672387565806504265/618970019642690137449562111 = 179622968672387565806504265 /( 2^89 - 1)
External Angle as a floating point decimal number = 2.9019655713870868535821260055542440298749779423213948304299730531995503353103626302473331181359966368582651105245850405837027542373052381532777325121338632071561064451614697645709384232759475708007812e-1
external angle as a binary rational (string) : 1001010010010100101001001010010010100101001001010010100100101001001010010100100101001001/11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
external angle as a binary floating number in exponential form =0.10010100100101001010010010100100101001010010010100101001001010010010100101001001010010010100101001001010010100100101001001010010100100101001010010010100100101001010010010100100101001010010010100101001*2^-1
external angle as a binary floating number in periodic form =0.(01001010010010100101001001010010010100101001001010010100100101001001010010100100101001001)
.(01001010010010100101001001010010010100101001001010010100100101001001010010100100101001001)
*/
#include <stdlib.h> // malloc
#include <stdio.h>
#include <gmp.h> // for rational numbers
#include <mpfr.h> // for floating point mumbers
// rotation map
//the number n is always increased by n0 modulo d
// input : op = n/d ( rational number ) and n0 ( integer)
// n = (n + n0 ) % d
// d = d
// output = rop = n/d
void mpq_rotation(mpq_t rop, const mpq_t op, const mpz_t n0)
{
mpz_t n; // numerator
mpz_t d; // denominator
mpz_inits( n, d, NULL);
//
mpq_get_num (n, op); //
mpq_get_den (d, op);
// n = (n + n0 ) % d
mpz_add(n, n, n0);
mpz_mod( n, n, d);
// output
mpq_set_num(rop, n);
mpq_set_den(rop, d);
mpz_clears( n, d, NULL);
}
void mpq_wake(mpq_t rop, mpq_t op)
{
// arbitrary precision variables from GMP library
mpz_t n0 ; // numerator of q
mpz_t nc;
mpz_t n;
mpz_t d ; // denominator of q
mpz_t m; // 2^i
mpz_t num ; // numerator of rop
mpz_t den ; // denominator of rop
long long int i;
unsigned long int base = 2;
unsigned long int id;
int cmp;
mpz_inits(n, n0,nc,d,num,den,m, NULL);
mpq_get_num(n0,op);
mpq_get_den(d,op);
id = mpz_get_ui(d);
// if (n <= 0 || n >= d ) error !!!! bad input
mpz_sub(nc, d, n0); // nc = d - n0
mpz_set(n, n0);
mpz_set_ui(num, 0);
// rop
// num = numerator(rop)
// denominator = den(rop) = (2^i) -1
mpz_ui_pow_ui(den, base, id) ; // den = base^id
mpz_sub_ui(den, den, 1); // den = den-1
// numerator
for (i=0; i<id ; i++){
mpz_set_ui(m, 0);
cmp = mpz_cmp(n,nc);// Compare op1 and op2. Return a positive value if op1 > op2, zero if op1 = op2, or a negative value if op1 < op2.
if ( cmp>0 ) {
mpz_ui_pow_ui(m, 2, id-i-1); // m = 2^(id-i )
mpz_add(num, num, m); // num = num + m
if (mpz_cmp(num, den) >0) mpz_mod( num, num, den); // num = num % d ; if num==d gives 0
//gmp_printf("s = 1");
}
// else gmp_printf("s = 0");
//gmp_printf (" i = %ld internal angle = %Zd / %Zd ea = %Zd / %Zd ; m = %Zd \n", i, n, d, num, den, m);
// n = (n + n0 ) % d = rotation
mpz_add(n, n, n0);
if (mpz_cmp(n, d)>0) mpz_mod( n, n, d);
//
//
}
// rop = external angle
mpq_set_num(rop,num);
mpq_set_den(rop,den);
mpq_canonicalize (rop); // It is the responsibility of the user to canonicalize the assigned variable before any arithmetic operations are performed on that variable.
// clear memory
mpz_clears(n, n0, nc, d, num,den, m, NULL);
}
/*
http://stackoverflow.com/questions/9895216/remove-character-from-string-in-c
"The idea is to keep a separate read and write pointers (pr for reading and pw for writing),
always advance the reading pointer, and advance the writing pointer only when it's not pointing to a given character."
modified
remove first length2rmv chars and after that take only length2stay chars from input string
output = input string
*/
void extract_str(char* str, unsigned int length2rmv, unsigned long int length2stay) {
// separate read and write pointers
char *pr = str; // read pointer
char *pw = str; // write pointer
int i =0; // index
while (*pr) {
if (i>length2rmv-1 && i <length2rmv+length2stay)
pw += 1; // advance the writing pointer only when
pr += 1; // always advance the reading pointer
*pw = *pr;
i +=1;
}
*pw = '\0';
}
int main ()
{
// notation :
//number type : s = string ; q = rational ; z = integer, f = floating point
// base : b = binary ; d = decimal
char *sqdInternalAngle = "13/34";
mpq_t qdInternalAngle; // internal angle = rational number q = n/d
mpz_t den;
unsigned long int uiIADenominator;
mpq_t qdExternalAngle; // rational number q = n/d
char *sqbExternalAngle;
mpfr_t fdExternalAngle ; //
char *sfbExternalAngle; //
mp_exp_t exponent ; // holds the exponent for the result string
mpz_t zdEANumerator;
mpz_t zdEADenominator;
mpfr_t EANumerator;
mpfr_t EADenominator;
mpfr_prec_t p = 200; // in bits , should be > denominator of internal angle
mpfr_set_default_prec (p); // but previously initialized variables are unaffected.
//mpfr_set_default_prec (precision);
// init variables
//mpf_init(fdExternalAngle);
mpz_inits(den, zdEANumerator,zdEADenominator, NULL);
mpq_inits (qdExternalAngle, qdInternalAngle, NULL); //
mpfr_inits(fdExternalAngle, EANumerator, EADenominator, NULL);
// set variables
mpq_set_str(qdInternalAngle, sqdInternalAngle, 10); // string is an internal angle
mpq_canonicalize (qdInternalAngle); // It is the responsibility of the user to canonicalize the assigned variable before any arithmetic operations are performed on that variable.
mpq_get_den(den,qdInternalAngle);
uiIADenominator = mpz_get_ui(den);
printf("uiIADenominator = %lu \n", uiIADenominator);
if ( p < uiIADenominator) printf("increase precision !!!!\n");
mpfr_printf("Using MPFR-%s with GMP-%s with precision = %u bits \n", mpfr_version, gmp_version, (unsigned int) p);
//
mpq_wake(qdExternalAngle, qdInternalAngle); // internal -> external
mpq_get_num(zdEANumerator ,qdExternalAngle);
mpq_get_den(zdEADenominator,qdExternalAngle);
// conversions
mpfr_set_z (EANumerator, zdEANumerator, GMP_RNDN);
mpfr_set_z (EADenominator, zdEADenominator, GMP_RNDN);
sqbExternalAngle = mpq_get_str (NULL, 2, qdExternalAngle); // rational number = fraction : from decimal to binary
mpfr_div (fdExternalAngle, EANumerator, EADenominator, GMP_RNDN);
sfbExternalAngle = (char*)malloc((sizeof(char) * uiIADenominator*2*4) + 3);
// mpfr_get_str (char *str, mpfr_exp_t *expptr, int b, size_t n, mpfr_t op, mpfr_rnd_t rnd)
if (sfbExternalAngle==NULL ) {printf("sfbExternalAngle error \n"); return 1;}
mpfr_get_str(sfbExternalAngle, &exponent, 2,200, fdExternalAngle, GMP_RNDN);
// print
gmp_printf ("internal angle = %Qd\n", qdInternalAngle); //
printf("first external angle : \n");
gmp_printf ("period = denominator of internal angle = %Zd\n", den); //
gmp_printf ("external angle as a decimal fraction = %Qd = %Zd /( 2^%Zd - 1) \n", qdExternalAngle, zdEANumerator, den); //
printf ("External Angle as a floating point decimal number = ");
mpfr_out_str (stdout, 10, p, fdExternalAngle, MPFR_RNDD); putchar ('\n');
gmp_printf ("external angle as a binary rational (string) : %s \n", sqbExternalAngle); //
printf ("external angle as a binary floating number in exponential form =0.%s*%d^%ld\n", sfbExternalAngle, 2, exponent);
extract_str(sfbExternalAngle, uiIADenominator+exponent, uiIADenominator);
printf ("external angle as a binary floating number in periodic form =0.(%s)\n", sfbExternalAngle);
// clear memory
//mpf_clear(fdExternalAngle);
mpq_clears(qdExternalAngle, qdInternalAngle, NULL);
mpz_clears(den, zdEANumerator, zdEADenominator, NULL);
mpfr_clears(fdExternalAngle, EANumerator, EADenominator, NULL);
free(sfbExternalAngle);
return 0;
}