Commit 168baa34 authored by SOLIMAN Sylvain's avatar SOLIMAN Sylvain

removed files that should not be here

parent 341327d1
CC=swipl-ld
LDFLAGS=-shared
SRC= $(wildcard *.c)
OBJ= $(SRC:.c=.o)
LDLIBS=-lgsl -lgslcblas -lm
$(foreach var, PLSOEXT, \
$(eval \
$(shell \
swipl -dump-runtime-variables | \
grep ^$(var)= | \
sed -E 's/="/=/;s/";$$//')))
all: roots test
.PHONY: clean
clean:
rm roots.o
roots: $(OBJ)
$(CC) $(LDFLAGS) -o $@ $^ $(LDLIBS)
mv $@.$(PLSOEXT) $@
roots.o: roots.h
%.o: %.c
$(CC) -o $@ -c $< $(CFLAGS)
test:
LD_LIBRARY_PATH=$$LD_LIBRARY_PATH$${LD_LIBRARY_PATH:+:}$$PWD
\ No newline at end of file
#include <stdio.h>
#include <stdlib.h>
#include <SWI-Prolog.h>
#include <gsl/gsl_poly.h>
#include "partialfraction.h"
static int
near(const double x, const double y)
{
return (x - y < EPSILON_POL) && (y - x < EPSILON_POL);
}
static void
system_solve(const double matrix_data[],
const double vector_data[],
const int n,
gsl_vector * x)
{
int s;
gsl_matrix_view m
= gsl_matrix_view_array(matrix_data, n, n);
gsl_vector_view b
= gsl_vector_view_array(vector_data, n);
gsl_permutation * p = gsl_permutation_alloc(n);
gsl_linalg_LU_decomp(&m.matrix, p, &s);
gsl_linalg_LU_solve(&m.matrix, p, &b.vector, x);
gsl_permutation_free(p);
}
/* factors a polynomial into a product of the form
(x + alpha_0)^n_0 ... (x2 + beta_0 x + gamma_0)^m_0 */
static void
get_factors(const double z[], const int N,
double ** alpha, double ** beta, double ** gamma,
int ** n, int ** m,
int * p, int * q)
{
int i, j, k;
int * cl, * re;
*p = *q = 0;
cl = malloc(N * sizeof(int));
re = malloc(N * sizeof(int));
for(i = 0 ; i < N ; i++)
cl[i] = 0;
for(i = 0 ; i < N ; i++)
{
/* new root encountered */
if(cl[i] == 0) {
/* real root */
if(near(z[2*i+1], 0))
{
cl[i] = ++(*p);
re[i] = 1;
}
/* imaginary root */
else
{
cl[i] = ++(*q);
re[i] = 0;
}
/* look up for copies of the root and its conjugate */
for(j = i + 1 ; j < N ; j++)
{
if( (cl[j] == 0) &&
near(z[2*i], z[2*j]) &&
( near(z[2*i+1], z[2*j+1]) || near(z[2*i+1], -z[2*j+1]) ) )
{
cl[j] = cl[i];
re[j] = re[i];
}
}
}
}
*alpha = malloc(*p * sizeof(double));
*beta = malloc(*q * sizeof(double));
*gamma = malloc(*q * sizeof(double));
*n = malloc(*p * sizeof(int));
*m = malloc(*q * sizeof(int));
/* count the multiplicity of the real roots +
set coefficient alpha */
for(k = 0 ; k < *p ; k++)
{
(*n)[k] = 0;
for(i = 0 ; i < N ; i++)
{
if(cl[i] == k + 1 && re[i] == 1)
{
if((*n)[k] == 0)
(*alpha)[k] = -z[2*i];
(*n)[k]++;
}
}
}
/* count the multiplicity of the imaginary roots +
set coefficients beta and gamma */
for(k = 0 ; k < *q ; k++)
{
(*m)[k] = 0;
for(i = 0 ; i < N ; i++)
{
if(cl[i] == k + 1 && re[i] == 0)
{
if((*m)[k] == 0)
{
(*beta)[k] = -2*z[2*i];
(*gamma)[k] = z[2*i]*z[2*i] + z[2*i+1]*z[2*i+1];
}
(*m)[k]++;
}
}
/* divide by two since every im. root is encountered twice */
(*m)[k] /= 2;
}
}
static foreign_t
pl_roots (term_t l,// list of coefficients
term_t s,// size of the list
term_t r // list of roots
)
{
int _n;
int i;
double * a, * z;
double * alpha, * beta, * gamma;
int p, q;
int * n, * m;
term_t head = PL_new_term_ref();
term_t list = PL_copy_term_ref(l);
term_t root = PL_new_term_ref();
if(!PL_get_integer(s, &_n))
PL_fail;
if(_n == 0)
PL_fail;
a = malloc(_n * sizeof(double));
z = malloc(2*(_n-1) * sizeof(double));
i = 0;
while(PL_get_list(list, head, list) && (i < _n))
{
double c;
if(PL_get_float(head, &c))
a[i] = c;
else
PL_fail;
i++;
}
gsl_poly_complex_workspace * w
= gsl_poly_complex_workspace_alloc (_n);
gsl_poly_complex_solve (a, _n, w, z);
gsl_poly_complex_workspace_free (w);
get_factors(z, _n-1, &alpha, &beta, &gamma, &n, &m, &p, &q);
free(a);
free(z);
free(alpha);
free(beta);
free(gamma);
free(n);
free(m);
return PL_unify_nil(r);
}
static foreign_t
pl_display(term_t l)
{
term_t head = PL_new_term_ref();
term_t list = PL_copy_term_ref(l);
while(PL_get_list(list, head, list))
{
double a;
if(PL_get_float(head, &a))
printf("%lf\n", a);
else
PL_fail;
}
return PL_get_nil(list);
}
install_t
install_partialfraction()
{
PL_register_foreign("roots", 3, pl_roots, 0);
PL_register_foreign("display", 1, pl_display, 0);
}
#ifndef PARTIALFRACTION_H
#define PARTIALFRACTION_H
#define EPSILON_POL 0.0001
#endif
:- module(partialfraction, [
partial_fraction/5
]).
:- use_foreign_library(foreign(roots)).
% A polynomial is represented by the list of its coefficients
% the fist one is the constant and the last one is the leading coefficient
% The null polynomial is represented by the empty list
% A polynomial is under its proper form if its l.c. is not zero.
% this predicate is used instead of the test to zero
small(A) :-
A < 0.0001, A > -0.0001.
% calculates the degree of a polynomial
% with the null polynomial being of degree -1 instead of -infinity.
degree([], -1).
degree([_|A_next], D) :-
degree(A_next, D_next),
D is D_next + 1.
% calculates the oppisite polynomial, i.e. -P.
opposite([], []).
opposite([A|A_next], [B|B_next]) :-
B is -A,
opposite(A_next, B_next).
% calculates the leading coefficient of a polynomial.
leading([], 0).
leading([A|A_next], LC) :-
leading(A_next, X),
( X =:= 0 -> LC is A ; LC is X ).
% calculates the difference betwenn two polynomial, possibly not in proper form.
row_subtract(A, [], A) :- !.
row_subtract([], B, C) :-
opposite(B, C).
row_subtract([A|A_next], [B|B_next], [C|C_next]) :-
C is float(A - B),
row_subtract(A_next, B_next, C_next).
% calculates the sum of two polynomials, possibly not in proper form.
row_add(A, [], A) :- !.
row_add([], B, B) :- !.
row_add([A|A_next], [B|B_next], [C|C_next]) :-
C is float(A + B),
row_add(A_next, B_next, C_next).
% multiply a polynomial by a constant C.
mult(0, _, []) :- !.
mult(_, [], []) :- !.
mult(C, [A0|A_next], [B0|B_next]) :-
B0 is C * A0,
mult(C, A_next, B_next).
% calculates the product of a polynomial by a given monomial
% whose degree is D and coefficient is C.
monomial_mult(0, C, A, B) :-
mult(C, A, B), !.
monomial_mult(D, C, A, [0|B]) :-
Dm is D - 1,
monomial_mult(Dm, C, A, B).
% delete zeros at the end, i.e. calculates the proper form
simplify([], []).
simplify([N|T], []) :-
small(N),
simplify(T, []), !.
simplify([A|Anxt], [A|Asmpl]) :-
simplify(Anxt, Asmpl).
% calculates the unitary polynomial associated
normalized([], []) :- !.
normalized(A, Anorm) :-
leading(A, Lc),
mult(float(1 rdiv Lc), A, Anorm).
% subtracts two polynomials
subtract(A, B, C) :-
row_subtract(A, B, D),
simplify(D, C).
% adds two polynomials
add(A, B, C) :-
row_add(A, B, D),
simplify(D, C).
% multiply two polynomials
prod(_, [], []) :- !.
prod(A, [B|Bnxt], C) :-
mult(B, A, C0),
prod(A, Bnxt, C1), !,
add(C0, [0|C1], C).
% calculates Q and R s.t. A = BQ + R with deg(R) < deg(B).
% [in floating point arithmetic]
long_division(_, [], _, _) :-
fail.
long_division(A, B, Q, R) :-
degree(A, DegA),
degree(B, DegB),
(
DegA < DegB
->
Q = [],
R = A
;
leading(A, LcA),
leading(B, LcB),
D is DegA - DegB,
C is float(LcA rdiv LcB),
monomial_mult(D, C, B, B1),
subtract(A, B1, A1),
long_division(A1, B, Q1, R), !,
monomial_mult(D, C, [1], Q0),
add(Q0, Q1, Q)
).
% calculates a greatest common divisor D of two polynomials
% and (U, V) s.t. AU + BV = D
euclid(A, [], [1], [], A) :- !.
euclid([], A, [], [1], A) :- !.
euclid(A, B, U, V, D) :-
degree(A, DegA),
degree(B, DegB),
DegA >= DegB, !,
long_division(A, B, Q, R),
euclid(B, R, Uprime, Vprime, D),
U = Vprime,
prod(Q, Vprime, T),
subtract(Uprime, T, V).
euclid(A, B, U, V, D) :-
euclid(B, A, V, U, D).
% reduces a rational fraction into a fraction whose
% denominator and numerator are coprimes.
simplify_fraction(A, B, Aprime, Bprime) :-
euclid(A, B, _, _, D),
long_division(A, D, Aprime, _),
long_division(B, D, Bprime, _).
% calculates the only unitary gcd.
normalized_euclid(A, B, Unorm, Vnorm, Dnorm) :-
euclid(A, B, U, V, D),
leading(D, LcD),
C is float(1 rdiv LcD),
mult(C, U, Unorm),
mult(C, V, Vnorm),
mult(C, D, Dnorm), !.
prod_list([], [1]).
prod_list([Q|L], P) :-
prod_list(L, R),
prod(Q, R, P).
exp(_, 0, [1]) :- !.
exp(P, N, Q) :-
Nprime is N - 1,
exp(P, Nprime, R),
prod(P, R, Q).
% builds a polynomial from a list of factors with multiplicities.
prod_factors([], [], [1]).
prod_factors([P|PT], [M|MT], Q) :-
prod_factors(PT, MT, R),
exp(P, M, S),
prod(R, S, Q).
% builds a polynomial from a list of factors with multiplicities
% except for the nth polynomial of the list.
prod_special([_|PT], [_|MT], 0, Q) :-
prod_factors(PT, MT, Q), !.
prod_special([P|PT], [M|MT], N, Q) :-
Nprime is N - 1,
prod_special(PT, MT, Nprime, R),
exp(P, M, S),
prod(R, S, Q).
% calculates the product of a list of polynomials by a given polynomial.
prod_list([], _, []).
prod_list([P|PT], Q, [R|RT]) :-
prod(P, Q, R),
prod_list(PT, Q, RT).
dotprod_list([], [], []).
dotprod_list([P|PT], [Q|QT], [R|RT]) :-
prod(P, Q, R),
dotprod_list(PT, QT, RT).
% gives bezout relation for a list of polynomials
general_bezout([A], [[1]], A) :- !.
general_bezout([P|PT], [U|QT], D) :-
general_bezout(PT, RT, E),
normalized_euclid(P, E, U, V, D),
prod_list(RT, V, QT).
% returns the list of the co-polynomials of a list of polynomial
% i.e. the list [Q/Q-0, Q/Q_1, ...] where Q = Q_0 ... Q_n.
co_polynomials(QL, ML, QLco) :-
built_co_polynomials(QL, ML, 0, QLco).
built_co_polynomials(QL, _, N, []) :-
length(QL, N), !.
built_co_polynomials(QL, ML, N, [Qco|QTco]) :-
prod_special(QL, ML, N, Qco),
Nnext is N + 1,
built_co_polynomials(QL, ML, Nnext, QTco).
% from P/Q_0...Q_n get S0/Q0^m_0 + ... + Sn/Qn^m_n
pre_partial_fraction(P, QL, ML, SL) :-
co_polynomials(QL, ML, TL),
general_bezout(TL, UL, _),
prod_list(UL, P, SL).
% from S/Q^m gets A_m/Q^m + ... + A_1/Q with deg(A_i) < deg(Q)
local_partial_fraction(_, _, 0, []) :- !.
local_partial_fraction(S, Q, M, [A|AT]) :-
long_division(S, Q, Sprime, A),
Mprime is M - 1,
local_partial_fraction(Sprime, Q, Mprime, AT).
list_local_partial_fraction([], [], [], []) :- !.
list_local_partial_fraction([S|ST], [Q|QT], [M|MT], [A|AT]) :-
local_partial_fraction(S, Q, M, A),
list_local_partial_fraction(ST, QT, MT, AT).
% calculates partial fraction decomposition from numerator
% and denominator decomposition.
partial_fraction(P, QL, ML, AL) :-
pre_partial_fraction(P, QL, ML, SL),
list_local_partial_fraction(SL, QL, ML, AL).
/* Inputs:
P, Q,
Outputs:
A is a list of polynomials of the form A_j_i/Q_i^j (i increasing, j decreasing),
QL is a list of irreductible unitary factors,
ML is the list of multiplicities with respect to QL.
*/
partial_fraction(P, Q, A, QL, ML) :-
roots(Q, QLraw, ML),
clean_recursive(QLraw, QL),
partial_fraction(P, QL, ML, Araw),
clean_recursive(Araw, A).
% to compensate rounding errors and stay clean
clean_number(X, Y) :-
Xr is round(X),
ERR is X - Xr,
( small(ERR) -> Y = Xr ; Y = X).
clean_recursive(X, Y) :-
(X = [] -> Y = []
;
X = [H|T] -> clean_recursive(H, CH),
clean_recursive(T, CT),
Y = [CH|CT]
;
clean_number(X, Y)
).
#include <stdio.h>
#include <stdlib.h>
#include <SWI-Prolog.h>
#include <gsl/gsl_poly.h>
#include "roots.h"
static int
near(const double x, const double y)
{
return (x - y < EPSILON_POL) && (y - x < EPSILON_POL);
}
/* static void
system_solve(const double matrix_data[],
const double vector_data[],
const int n,
gsl_vector * x)
{
int s;
gsl_matrix_view m
= gsl_matrix_view_array(matrix_data, n, n);
gsl_vector_view b
= gsl_vector_view_array(vector_data, n);
gsl_permutation * p = gsl_permutation_alloc(n);
gsl_linalg_LU_decomp(&m.matrix, p, &s);
gsl_linalg_LU_solve(&m.matrix, p, &b.vector, x);
gsl_permutation_free(p);
}*/
/* factors a polynomial into a product of the form
(x + alpha_0)^n_0 ... (x2 + beta_0 x + gamma_0)^m_0 */
static void
get_factors(const double z[], const int N,
double ** alpha, double ** beta, double ** gamma,
int ** n, int ** m,
int * p, int * q)
{
int i, j, k;
int * cl, * re;
*p = *q = 0;
cl = malloc(N * sizeof(int));
re = malloc(N * sizeof(int));
for(i = 0 ; i < N ; i++)
cl[i] = 0;
for(i = 0 ; i < N ; i++)
{
/* new root encountered */
if(cl[i] == 0) {
/* real root */
if(near(z[2*i+1], 0))
{
cl[i] = ++(*p);
re[i] = 1;
}
/* imaginary root */
else
{
cl[i] = ++(*q);
re[i] = 0;
}
/* look up for copies of the root and its conjugate */
for(j = i + 1 ; j < N ; j++)
{
if( (cl[j] == 0) &&
near(z[2*i], z[2*j]) &&
( near(z[2*i+1], z[2*j+1]) || near(z[2*i+1], -z[2*j+1]) ) )
{
cl[j] = cl[i];
re[j] = re[i];
}
}
}
}
*alpha = malloc(*p * sizeof(double));
*beta = malloc(*q * sizeof(double));
*gamma = malloc(*q * sizeof(double));
*n = malloc(*p * sizeof(int));
*m = malloc(*q * sizeof(int));
/* count the multiplicity of the real roots +
set coefficient alpha */
for(k = 0 ; k < *p ; k++)
{
(*n)[k] = 0;
for(i = 0 ; i < N ; i++)
{
if(cl[i] == k + 1 && re[i] == 1)
{
if((*n)[k] == 0)
(*alpha)[k] = -z[2*i];
(*n)[k]++;
}
}
}
/* count the multiplicity of the imaginary roots +
set coefficients beta and gamma */
for(k = 0 ; k < *q ; k++)
{
(*m)[k] = 0;
for(i = 0 ; i < N ; i++)
{
if(cl[i] == k + 1 && re[i] == 0)
{
if((*m)[k] == 0)
{
(*beta)[k] = -2*z[2*i];
(*gamma)[k] = z[2*i]*z[2*i] + z[2*i+1]*z[2*i+1];
}
(*m)[k]++;
}
}
/* divide by two since every im. root is encountered twice */
(*m)[k] /= 2;
}
}
/* associate a polynomial with the list of its irreducible factors */
static foreign_t
pl_roots (term_t l,// list of coefficients
term_t irr, // list of irreducible polynomials
term_t mul // list of multiplicities
)
{
int _n = 0;
int i = 0, j;
double * a, * z;
double * alpha, * beta, * gamma;
int p, q;
int * n, * m;
term_t head = PL_new_term_ref();
term_t list = PL_copy_term_ref(l);
term_t factor = PL_new_term_ref();
term_t coeff = PL_new_term_ref();
term_t expst = PL_new_term_ref();
/* read the list of coefficients */
while(PL_get_list(list, head, list))
_n++;
list = PL_copy_term_ref(l);
a = malloc(_n * sizeof(double));
z = malloc(2 * (_n - 1) * sizeof(double));
while(PL_get_list(list, head, list))
{
double c;
if(PL_get_float(head, &c))
a[i] = c;
else
PL_fail;
i++;
}
/* pathological cases */
if(_n < 2)
PL_fail;