![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to calcul_integral.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / src|
Return to calcul_integral.c CVS log | Up to [local] / rpl / src |
1.1 bertrand 1: /*
2: ================================================================================
1.22 ! bertrand 3: RPL/2 (R) version 4.1.0.prerelease.2
1.15 bertrand 4: Copyright (C) 1989-2011 Dr. BERTRAND Joël
1.1 bertrand 5:
6: This file is part of RPL/2.
7:
8: RPL/2 is free software; you can redistribute it and/or modify it
9: under the terms of the CeCILL V2 License as published by the french
10: CEA, CNRS and INRIA.
11:
12: RPL/2 is distributed in the hope that it will be useful, but WITHOUT
13: ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14: FITNESS FOR A PARTICULAR PURPOSE. See the CeCILL V2 License
15: for more details.
16:
17: You should have received a copy of the CeCILL License
18: along with RPL/2. If not, write to info@cecill.info.
19: ================================================================================
20: */
21:
22:
1.11 bertrand 23: #include "rpl-conv.h"
1.1 bertrand 24:
25:
26: /*
27: ================================================================================
28: Fonction 'integrale_romberg'
29: ================================================================================
30: Entrées : pointeur sur une structure struct_processus
31: --------------------------------------------------------------------------------
32: Sorties :
33: --------------------------------------------------------------------------------
34: Effets de bord : néant
35: ================================================================================
36: */
37:
38: void
39: integrale_romberg(struct_processus *s_etat_processus,
40: struct_objet *s_expression, unsigned char *variable,
41: real8 a, real8 b, real8 precision)
42: {
43: real8 erreur;
44: real8 h;
45: real8 **t;
46: real8 **t_tampon;
47: real8 *vecteur;
48: real8 *vecteur_precedent;
49: real8 x;
50:
51: logical1 validite;
52: logical1 erreur_memoire;
53:
54: struct_objet *s_objet;
55:
56: struct_variable s_variable;
57:
58: unsigned long i;
59: unsigned long j;
60: unsigned long k;
61: unsigned long n;
62: unsigned long p;
63: unsigned long nombre_elements;
64: unsigned long nombre_termes;
65: unsigned long taille_vecteur_precedent;
66:
67: /*
68: * Création d'une variable locale représentant la variable d'intégration
69: */
70:
71: (*s_etat_processus).niveau_courant++;
72: s_variable.niveau = (*s_etat_processus).niveau_courant;
73:
74: if ((s_variable.nom = (unsigned char *) malloc((strlen(variable) + 1) *
75: sizeof(unsigned char))) == NULL)
76: {
77: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
78: return;
79: }
80:
81: strcpy(s_variable.nom, variable);
82:
83: if ((s_variable.objet = allocation(s_etat_processus, REL)) == NULL)
84: {
85: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
86: return;
87: }
88:
89: if (creation_variable(s_etat_processus, &s_variable, 'V', 'P') == d_erreur)
90: {
91: return;
92: }
93:
94: /*
95: * Algorithme de Romberg
96: */
97:
98: if ((vecteur = (real8 *) malloc(2 * sizeof(real8))) == NULL)
99: {
100: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
101: return;
102: }
103:
104: evaluation_romberg(s_etat_processus, s_expression, variable, &a,
105: &(vecteur[nombre_termes = 0]), &validite);
106: nombre_termes += (validite == d_vrai) ? 1 : 0;
107:
108: evaluation_romberg(s_etat_processus, s_expression, variable, &b,
109: &(vecteur[nombre_termes]), &validite);
110: nombre_termes += (validite == d_vrai) ? 1 : 0;
111:
112: if ((t = (real8 **) malloc(sizeof(real8 *))) == NULL)
113: {
114: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
115: return;
116: }
117:
118: if ((t[0] = (real8 *) malloc(sizeof(real8))) == NULL)
119: {
120: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
121: return;
122: }
123:
124: t[0][n = 0] = ((b - a) / 2) * sommation_vecteur_reel(vecteur,
125: &nombre_termes, &erreur_memoire);
126:
127: if (erreur_memoire == d_vrai)
128: {
129: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
130: return;
131: }
132:
133: free(vecteur);
134:
135: vecteur_precedent = NULL;
136: taille_vecteur_precedent = 0;
137:
138: do
139: {
140: h = (b - a) / (nombre_elements = (1ULL << (++n)));
141: x = a;
142:
143: /*
144: * Nouvelle allocation de t
145: */
146:
147: t_tampon = t;
148:
149: if ((t = (real8 **) malloc((n + 1) * sizeof(real8 *))) == NULL)
150: {
151: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
152: return;
153: }
154:
155: for(i = 0; i <= n; i++)
156: {
157: if ((t[i] = (real8 *) malloc((n + 1) * sizeof(real8))) == NULL)
158: {
159: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
160: return;
161: }
162:
163: if (i != n)
164: {
165: for(t[i][n] = 0, j = 0; j < n; j++)
166: {
167: t[i][j] = t_tampon[i][j];
168: }
169: }
170: else
171: {
172: for(j = 0; j <= n; t[i][j++] = 0);
173: }
174: }
175:
176: for(i = 0; i < n; free(t_tampon[i++]));
177: free(t_tampon);
178:
179: /*
180: * Boucle principale
181: */
182:
183: if ((vecteur = (real8 *) malloc((nombre_elements + 1) *
184: sizeof(real8))) == NULL)
185: {
186: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
187: return;
188: }
189:
190:
191: if (vecteur_precedent == NULL)
192: {
193: /*
194: * Première itération
195: */
196:
197: evaluation_romberg(s_etat_processus, s_expression, variable, &a,
198: &(vecteur[nombre_termes = 0]), &validite);
199: nombre_termes += (validite == d_vrai) ? 1 : 0;
200:
201: evaluation_romberg(s_etat_processus, s_expression, variable, &b,
202: &(vecteur[nombre_termes]), &validite);
203: nombre_termes += (validite == d_vrai) ? 1 : 0;
204:
205: if (nombre_termes == 2)
206: {
207: vecteur[0] += vecteur[1];
208: vecteur[0] /= 2;
209: nombre_termes--;
210: }
211:
212: for(i = 1; i < nombre_elements; i++)
213: {
214: x += h;
215:
216: evaluation_romberg(s_etat_processus, s_expression, variable, &x,
217: &(vecteur[nombre_termes]), &validite);
218: nombre_termes += (validite == d_vrai) ? 1 : 0;
219: }
220: }
221: else
222: {
223: /*
224: * Pour les itérations suivantes, la moitié des points
225: * a déjà été calculée.
226: */
227:
228: for(i = 0; i < taille_vecteur_precedent; i++)
229: {
230: vecteur[i] = vecteur_precedent[i];
231: }
232:
233: free(vecteur_precedent);
234:
235: for(nombre_termes = taille_vecteur_precedent, i = 1;
236: i < nombre_elements; x += h, i += 2)
237: {
238: x += h;
239:
240: evaluation_romberg(s_etat_processus, s_expression, variable, &x,
241: &(vecteur[nombre_termes]), &validite);
242: nombre_termes += (validite == d_vrai) ? 1 : 0;
243: }
244: }
245:
246: t[0][n] = h * sommation_vecteur_reel(vecteur, &nombre_termes,
247: &erreur_memoire);
248:
249: if (erreur_memoire == d_vrai)
250: {
251: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
252: return;
253: }
254:
255: vecteur_precedent = vecteur;
256: taille_vecteur_precedent = nombre_termes;
257:
258: for(p = 1; p <= n; p++)
259: {
260: k = n - p;
261: t[p][k] = ((pow(4, p) * t[p - 1][k + 1]) - t[p - 1][k]) /
262: ((real8) (pow(4, p) - 1));
263: }
264: } while(((erreur = fabs(t[n][0] - t[n - 1][0])) > precision) &&
265: ((*s_etat_processus).var_volatile_requete_arret == 0));
266:
267: free(vecteur_precedent);
268:
269: /*
270: * Empilement du résultat et liberation de t
271: */
272:
273: if ((s_objet = allocation(s_etat_processus, REL)) == NULL)
274: {
275: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
276: return;
277: }
278:
279: (*((real8 *) (*s_objet).objet)) = t[n][0];
280:
281: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
282: s_objet) == d_erreur)
283: {
284: return;
285: }
286:
287: if ((s_objet = allocation(s_etat_processus, REL)) == NULL)
288: {
289: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
290: return;
291: }
292:
293: (*((real8 *) (*s_objet).objet)) = erreur;
294:
295: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
296: s_objet) == d_erreur)
297: {
298: return;
299: }
300:
301: for(i = 0; i <= n; free(t[i++]));
302: free(t);
303:
304: /*
305: * Libération de la variable locale
306: */
307:
308: (*s_etat_processus).niveau_courant--;
309:
310: if (retrait_variable(s_etat_processus, s_variable.nom, 'L')
311: == d_erreur)
312: {
313: if ((*s_etat_processus).erreur_systeme != d_es)
314: {
315: if ((*s_etat_processus).erreur_systeme ==
316: d_es_variable_introuvable)
317: {
318: (*s_etat_processus).erreur_systeme = d_es;
319: }
320: else
321: {
322: /*
323: * Erreur système
324: */
325:
326: return;
327: }
328: }
329:
330: (*s_etat_processus).erreur_execution = d_ex_variable_non_definie;
331: return;
332: }
333:
334: return;
335: }
336:
337:
338: void
339: evaluation_romberg(struct_processus *s_etat_processus,
340: struct_objet *s_expression, unsigned char *variable, real8 *point,
341: real8 *valeur, logical1 *validite)
342: {
343: long hauteur_pile;
344:
345: struct_liste_pile_systeme *l_position_normale;
346:
347: struct_objet *s_objet;
348:
349: /*
350: * Vérification du type de la variable
351: */
352:
353: (*validite) = d_faux;
354: (*valeur) = 0;
355:
356: if (recherche_variable(s_etat_processus, variable) == d_vrai)
357: {
1.18 bertrand 358: if ((*(*s_etat_processus).pointeur_variable_courante)
359: .variable_verrouillee == d_vrai)
1.1 bertrand 360: {
361: (*s_etat_processus).erreur_execution = d_ex_variable_verrouillee;
362: return;
363: }
364:
1.18 bertrand 365: if ((*(*(*s_etat_processus).pointeur_variable_courante)
1.1 bertrand 366: .objet).type != REL)
367: {
1.18 bertrand 368: liberation(s_etat_processus, (*(*s_etat_processus)
369: .pointeur_variable_courante).objet);
1.1 bertrand 370:
1.18 bertrand 371: if (((*(*s_etat_processus).pointeur_variable_courante).objet =
1.1 bertrand 372: allocation(s_etat_processus, REL)) == NULL)
373: {
374: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
375: return;
376: }
377: }
378:
1.18 bertrand 379: (*((real8 *) (*(*(*s_etat_processus).pointeur_variable_courante)
1.1 bertrand 380: .objet).objet)) = (*point);
381: }
382: else
383: {
384: /*
385: * La variable d'intégration étant locale, l'utilisateur ne peut
386: * pas l'effacer. On ne doit donc jamais passer par ici. Si
387: * c'est néanmoins le cas, une erreur système est générée
388: * provoquant l'arrêt du programme même dans une
389: * structure IFERR.
390: */
391:
392: return;
393: }
394:
395: hauteur_pile = (*s_etat_processus).hauteur_pile_operationnelle;
396: l_position_normale = (*s_etat_processus).l_base_pile_systeme;
397:
398: (*s_etat_processus).erreur_execution = d_ex;
399: (*s_etat_processus).exception = d_ep;
400:
401: if (evaluation(s_etat_processus, s_expression, 'N') == d_erreur)
402: {
403: if ((*s_etat_processus).erreur_systeme != d_es)
404: {
405: /*
406: * Erreur système
407: */
408:
409: return;
410: }
411: else
412: {
413: /*
414: * Restauration de la pile initiale
415: */
416:
417: while((*s_etat_processus).hauteur_pile_operationnelle >
418: (unsigned long) hauteur_pile)
419: {
420: if (depilement(s_etat_processus, &((*s_etat_processus)
421: .l_base_pile), &s_objet) == d_erreur)
422: {
423: (*s_etat_processus).erreur_execution =
424: d_ex_manque_argument;
425: break;
426: }
427:
428: liberation(s_etat_processus, s_objet);
429: }
430:
431: /*
432: * Restauration de la pile système
433: */
434:
435: while((*s_etat_processus).l_base_pile_systeme !=
436: l_position_normale)
437: {
438: depilement_pile_systeme(s_etat_processus);
439:
440: if ((*s_etat_processus).erreur_systeme != d_es)
441: {
442: /*
443: * Une pile vide provoque une erreur système
444: */
445:
446: return;
447: }
448: }
449: }
450:
451: (*s_etat_processus).erreur_execution = d_ex;
452: (*s_etat_processus).exception = d_ep;
453: }
454: else
455: {
456: /*
457: * Retrait de la valeur de la pile
458: */
459:
460: if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
461: &s_objet) == d_erreur)
462: {
463: (*s_etat_processus).erreur_execution = d_ex_manque_argument;
464: return;
465: }
466:
467: if ((*s_objet).type == INT)
468: {
469: (*valeur) = (*((integer8 *) (*s_objet).objet));
470: (*validite) = d_vrai;
471: }
472: else if ((*s_objet).type == REL)
473: {
474: (*valeur) = (*((real8 *) (*s_objet).objet));
475: (*validite) = d_vrai;
476: }
477:
478: liberation(s_etat_processus, s_objet);
479: }
480:
481: return;
482: }
483:
484: // vim: ts=4