← Back to branch summary

~ubuntu-branches/ubuntu/raring/scilab/raring-proposed

« back to all changes in this revision

Viewing changes to modules/special_functions/sci_gateway/c/sci_legendre.c

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2012-08-30 14:42:38 UTC
  • mfrom: (1.4.7)
  • Revision ID: package-import@ubuntu.com-20120830144238-c1y2og7dbm7m9nig
Tags: 5.4.0-beta-3-1~exp1
* New upstream release
* Update the scirenderer dep
* Get ride of libjhdf5-java dependency

Show diffs side-by-side

added added

removed removed

Lines of Context:
modules/special_functions/sci_gateway/c/sci_legendre.c
16
16
#include "localization.h"
17
17
#include "Scierror.h"
18
18
/*--------------------------------------------------------------------------*/
19
 
extern void C2F(dxlegf)(double *dnu1, int *nudiff, int *mu1, int *mu2, 
20
 
                        double *x,int *id, double *pqa, int *ipqa, int *ierror);
 
19
extern void C2F(dxlegf)(double *dnu1, int *nudiff, int *mu1, int *mu2,
 
20
                        double *x, int *id, double *pqa, int *ipqa, int *ierror);
21
21
/*--------------------------------------------------------------------------*/
22
22
static double return_an_inf(void);
23
23
static int verify_cstr(double x[], int nb_elt, int *xmin, int *xmax);
24
24
/*--------------------------------------------------------------------------*/
25
 
int sci_legendre(char *fname,unsigned long fname_len)
 
25
int sci_legendre(char *fname, unsigned long fname_len)
26
26
{
27
27
    /*
28
28
    *   Interface onto the (Slatec) dxleg.f code.
50
50
    double *x = NULL, xx = 0., dnu1 = 0., *pqa = NULL;
51
51
    int id = 0, ierror = 0, i = 0, j = 0, nudiff = 0;
52
52
 
53
 
    CheckLhs(1, 1); 
 
53
    CheckLhs(1, 1);
54
54
    CheckRhs(3, 4);
55
55
    GetRhsVar(1, MATRIX_OF_DOUBLE_DATATYPE, &mN, &nN, &lN);
56
56
 
57
 
    if ( ! verify_cstr(stk(lN), mN*nN, &n1, &n2) )
 
57
    if ( ! verify_cstr(stk(lN), mN * nN, &n1, &n2) )
58
58
    {
59
 
        Scierror(999,_("%s: Wrong type for first input argument.\n"), fname);
 
59
        Scierror(999, _("%s: Wrong type for first input argument.\n"), fname);
60
60
        return 0;
61
61
    };
62
62
 
63
63
    if ( mN == 1 && nN == 1) N_is_scalar = 1;
64
64
 
65
 
    GetRhsVar(2,MATRIX_OF_DOUBLE_DATATYPE, &mM, &nM, &lM);
66
 
    if ( ! verify_cstr(stk(lM), mM*nM, &m1, &m2) )
 
65
    GetRhsVar(2, MATRIX_OF_DOUBLE_DATATYPE, &mM, &nM, &lM);
 
66
    if ( ! verify_cstr(stk(lM), mM * nM, &m1, &m2) )
67
67
    {
68
 
        Scierror(999,_("%s: Wrong type for input argument #%d.\n"), fname,2);
 
68
        Scierror(999, _("%s: Wrong type for input argument #%d.\n"), fname, 2);
69
69
        return 0;
70
70
    }
71
71
 
73
73
 
74
74
    if ( ! M_is_scalar  &&  ! N_is_scalar )
75
75
    {
76
 
        Scierror(999,_("%s: Only one of arg1 and arg2 may be a vector.\n"), fname);
 
76
        Scierror(999, _("%s: Only one of arg1 and arg2 may be a vector.\n"), fname);
77
77
        return 0;
78
78
    };
79
79
 
80
 
    GetRhsCVar(3,MATRIX_OF_DOUBLE_DATATYPE, &it, &mx, &nx, &lx, &lc);
 
80
    GetRhsCVar(3, MATRIX_OF_DOUBLE_DATATYPE, &it, &mx, &nx, &lx, &lc);
81
81
    if ( it != 0 )
82
82
    {
83
 
        Scierror(999,_("%s: Wrong type for input argument #%d: Real matrix expected.\n"), fname, 3);
 
83
        Scierror(999, _("%s: Wrong type for input argument #%d: Real matrix expected.\n"), fname, 3);
84
84
        return 0;
85
85
    };
86
86
 
87
 
    mnx = mx*nx;
 
87
    mnx = mx * nx;
88
88
    x = stk(lx);
89
89
    for ( i = 0 ; i < mnx ; i++ )
90
90
        if ( ! (fabs(x[i]) < 1.0) )
91
91
        {
92
 
            Scierror(999,_("%s: Wrong value for input argument #%d: Matrix with elements in (%d,%d) expected.\n"), fname,3,-1,1);
 
92
            Scierror(999, _("%s: Wrong value for input argument #%d: Matrix with elements in (%d,%d) expected.\n"), fname, 3, -1, 1);
93
93
            return 0;
94
94
        };
95
95
 
96
96
    if ( Rhs == 4 )
97
97
    {
98
 
        GetRhsVar(4,STRING_DATATYPE, &ms, &ns, &ls);
99
 
        if ( strcmp(cstk(ls),"norm") == 0)
 
98
        GetRhsVar(4, STRING_DATATYPE, &ms, &ns, &ls);
 
99
        if ( strcmp(cstk(ls), "norm") == 0)
100
100
        {
101
101
            normalised = 1;
102
102
        }
112
112
 
113
113
    MNp1 = Max (n2 - n1, m2 - m1) + 1;
114
114
 
115
 
    CreateVar(Rhs+1, MATRIX_OF_DOUBLE_DATATYPE, &MNp1, &mnx, &lpqa); pqa = stk(lpqa);
116
 
    CreateVar(Rhs+2, MATRIX_OF_INTEGER_DATATYPE, &MNp1, &mnx, &lipqa); ipqa = istk(lipqa);
 
115
    CreateVar(Rhs + 1, MATRIX_OF_DOUBLE_DATATYPE, &MNp1, &mnx, &lpqa);
 
116
    pqa = stk(lpqa);
 
117
    CreateVar(Rhs + 2, MATRIX_OF_INTEGER_DATATYPE, &MNp1, &mnx, &lipqa);
 
118
    ipqa = istk(lipqa);
117
119
 
118
120
    if ( normalised )
119
121
    {
131
133
    {
132
134
        xx = fabs(x[i]); /* dxleg computes only for x in [0,1) */
133
135
        F2C(dxlegf) (&dnu1, &nudiff, &m1, &m2, &xx, &id,
134
 
            stk(lpqa+i*MNp1), istk(lipqa+i*MNp1), &ierror);
 
136
                     stk(lpqa + i * MNp1), istk(lipqa + i * MNp1), &ierror);
135
137
        if ( ierror != 0 )
136
138
        {
137
139
            if ( ierror == 207 ) /* @TODO what is 207 ? */
138
140
            {
139
 
                Scierror(999,_("%s: overflow or underflow of an extended range number\n"), fname);
 
141
                Scierror(999, _("%s: overflow or underflow of an extended range number\n"), fname);
140
142
            }
141
143
            else
142
144
            {
143
 
                Scierror(999,_("%s: error number %d\n"), fname, ierror);
 
145
                Scierror(999, _("%s: error number %d\n"), fname, ierror);
144
146
            }
145
147
            return 0;
146
148
        };
151
153
    *  When the "exponent" part (ipqa) is 0 then the number is exactly
152
154
    *  given by pqa else it leads to an overflow or an underflow.
153
155
    */
154
 
    for ( i = 0 ; i < mnx*MNp1 ; i++ )
 
156
    for ( i = 0 ; i < mnx * MNp1 ; i++ )
155
157
    {
156
158
        if ( ipqa[i] < 0 )
157
159
        {
164
166
    }
165
167
 
166
168
    /* complete the result by odd/even symmetry for negative x */
167
 
    for ( i = 0 ; i < mnx ; i++ ) 
 
169
    for ( i = 0 ; i < mnx ; i++ )
168
170
    {
169
 
        if ( x[i] < 0.0 ) 
 
171
        if ( x[i] < 0.0 )
170
172
        {
171
 
            if ( (n1+m1) % 2 == 1 ) 
 
173
            if ( (n1 + m1) % 2 == 1 )
172
174
            {
173
 
                for ( j = 0 ; j < MNp1 ; j+=2 )
 
175
                for ( j = 0 ; j < MNp1 ; j += 2 )
174
176
                {
175
 
                    pqa[i*MNp1 + j] = -pqa[i*MNp1 + j];
 
177
                    pqa[i * MNp1 + j] = -pqa[i * MNp1 + j];
176
178
                }
177
179
            }
178
 
            else 
 
180
            else
179
181
            {
180
 
                for ( j = 1 ; j < MNp1 ; j+=2 )
 
182
                for ( j = 1 ; j < MNp1 ; j += 2 )
181
183
                {
182
 
                    pqa[i*MNp1 + j] = -pqa[i*MNp1 + j];
 
184
                    pqa[i * MNp1 + j] = -pqa[i * MNp1 + j];
183
185
                }
184
186
            }
185
187
        }
199
201
 
200
202
    if ( first )
201
203
    {
202
 
        inf = inf/(inf - (double) first);
 
204
        inf = inf / (inf - (double) first);
203
205
        first = 0;
204
206
    }
205
207
    return (inf);
222
224
    }
223
225
    for ( i = 1 ; i < nb_elt ; i++ )
224
226
    {
225
 
        if ( x[i] != x[i-1]+1.0 )
 
227
        if ( x[i] != x[i - 1] + 1.0 )
226
228
        {
227
229
            return 0;
228
230
        }
229
231
    }
230
232
 
231
233
    *xmin = (int) x[0];
232
 
    *xmax = (int) x[nb_elt-1];
 
234
    *xmax = (int) x[nb_elt - 1];
233
235
    return 1;
234
236
}
235
237
/*--------------------------------------------------------------------------*/