![[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.42 ! bertrand 3: RPL/2 (R) version 4.1.13
1.41 bertrand 4: Copyright (C) 1989-2013 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)));
1.34 bertrand 141:
142: if (((unsigned long) nombre_elements) == 0)
143: {
144: // Dépassement de capacité
145: n--;
146: break;
147: }
148:
1.1 bertrand 149: x = a;
150:
151: /*
152: * Nouvelle allocation de t
153: */
154:
155: t_tampon = t;
156:
157: if ((t = (real8 **) malloc((n + 1) * sizeof(real8 *))) == NULL)
158: {
159: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
160: return;
161: }
162:
163: for(i = 0; i <= n; i++)
164: {
165: if ((t[i] = (real8 *) malloc((n + 1) * sizeof(real8))) == NULL)
166: {
167: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
168: return;
169: }
170:
171: if (i != n)
172: {
173: for(t[i][n] = 0, j = 0; j < n; j++)
174: {
175: t[i][j] = t_tampon[i][j];
176: }
177: }
178: else
179: {
180: for(j = 0; j <= n; t[i][j++] = 0);
181: }
182: }
183:
184: for(i = 0; i < n; free(t_tampon[i++]));
185: free(t_tampon);
186:
187: /*
188: * Boucle principale
189: */
190:
191: if ((vecteur = (real8 *) malloc((nombre_elements + 1) *
192: sizeof(real8))) == NULL)
193: {
194: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
195: return;
196: }
197:
198:
199: if (vecteur_precedent == NULL)
200: {
201: /*
202: * Première itération
203: */
204:
205: evaluation_romberg(s_etat_processus, s_expression, variable, &a,
206: &(vecteur[nombre_termes = 0]), &validite);
207: nombre_termes += (validite == d_vrai) ? 1 : 0;
208:
209: evaluation_romberg(s_etat_processus, s_expression, variable, &b,
210: &(vecteur[nombre_termes]), &validite);
211: nombre_termes += (validite == d_vrai) ? 1 : 0;
212:
213: if (nombre_termes == 2)
214: {
215: vecteur[0] += vecteur[1];
216: vecteur[0] /= 2;
217: nombre_termes--;
218: }
219:
220: for(i = 1; i < nombre_elements; i++)
221: {
222: x += h;
223:
224: evaluation_romberg(s_etat_processus, s_expression, variable, &x,
225: &(vecteur[nombre_termes]), &validite);
226: nombre_termes += (validite == d_vrai) ? 1 : 0;
227: }
228: }
229: else
230: {
231: /*
232: * Pour les itérations suivantes, la moitié des points
233: * a déjà été calculée.
234: */
235:
236: for(i = 0; i < taille_vecteur_precedent; i++)
237: {
238: vecteur[i] = vecteur_precedent[i];
239: }
240:
241: free(vecteur_precedent);
242:
243: for(nombre_termes = taille_vecteur_precedent, i = 1;
244: i < nombre_elements; x += h, i += 2)
245: {
246: x += h;
247:
248: evaluation_romberg(s_etat_processus, s_expression, variable, &x,
249: &(vecteur[nombre_termes]), &validite);
250: nombre_termes += (validite == d_vrai) ? 1 : 0;
251: }
252: }
253:
254: t[0][n] = h * sommation_vecteur_reel(vecteur, &nombre_termes,
255: &erreur_memoire);
256:
257: if (erreur_memoire == d_vrai)
258: {
259: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
260: return;
261: }
262:
263: vecteur_precedent = vecteur;
264: taille_vecteur_precedent = nombre_termes;
265:
266: for(p = 1; p <= n; p++)
267: {
268: k = n - p;
269: t[p][k] = ((pow(4, p) * t[p - 1][k + 1]) - t[p - 1][k]) /
270: ((real8) (pow(4, p) - 1));
271: }
272: } while(((erreur = fabs(t[n][0] - t[n - 1][0])) > precision) &&
273: ((*s_etat_processus).var_volatile_requete_arret == 0));
274:
275: free(vecteur_precedent);
276:
277: /*
278: * Empilement du résultat et liberation de t
279: */
280:
281: if ((s_objet = allocation(s_etat_processus, REL)) == NULL)
282: {
283: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
284: return;
285: }
286:
287: (*((real8 *) (*s_objet).objet)) = t[n][0];
288:
289: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
290: s_objet) == d_erreur)
291: {
292: return;
293: }
294:
295: if ((s_objet = allocation(s_etat_processus, REL)) == NULL)
296: {
297: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
298: return;
299: }
300:
301: (*((real8 *) (*s_objet).objet)) = erreur;
302:
303: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
304: s_objet) == d_erreur)
305: {
306: return;
307: }
308:
309: for(i = 0; i <= n; free(t[i++]));
310: free(t);
311:
312: /*
313: * Libération de la variable locale
314: */
315:
316: (*s_etat_processus).niveau_courant--;
317:
318: if (retrait_variable(s_etat_processus, s_variable.nom, 'L')
319: == d_erreur)
320: {
321: if ((*s_etat_processus).erreur_systeme != d_es)
322: {
323: if ((*s_etat_processus).erreur_systeme ==
324: d_es_variable_introuvable)
325: {
326: (*s_etat_processus).erreur_systeme = d_es;
327: }
328: else
329: {
330: /*
331: * Erreur système
332: */
333:
334: return;
335: }
336: }
337:
338: (*s_etat_processus).erreur_execution = d_ex_variable_non_definie;
339: return;
340: }
341:
342: return;
343: }
344:
345:
346: void
347: evaluation_romberg(struct_processus *s_etat_processus,
348: struct_objet *s_expression, unsigned char *variable, real8 *point,
349: real8 *valeur, logical1 *validite)
350: {
351: long hauteur_pile;
352:
353: struct_liste_pile_systeme *l_position_normale;
354:
355: struct_objet *s_objet;
356:
357: /*
358: * Vérification du type de la variable
359: */
360:
361: (*validite) = d_faux;
362: (*valeur) = 0;
363:
364: if (recherche_variable(s_etat_processus, variable) == d_vrai)
365: {
1.18 bertrand 366: if ((*(*s_etat_processus).pointeur_variable_courante)
367: .variable_verrouillee == d_vrai)
1.1 bertrand 368: {
369: (*s_etat_processus).erreur_execution = d_ex_variable_verrouillee;
370: return;
371: }
372:
1.18 bertrand 373: if ((*(*(*s_etat_processus).pointeur_variable_courante)
1.1 bertrand 374: .objet).type != REL)
375: {
1.18 bertrand 376: liberation(s_etat_processus, (*(*s_etat_processus)
377: .pointeur_variable_courante).objet);
1.1 bertrand 378:
1.18 bertrand 379: if (((*(*s_etat_processus).pointeur_variable_courante).objet =
1.1 bertrand 380: allocation(s_etat_processus, REL)) == NULL)
381: {
382: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
383: return;
384: }
385: }
386:
1.18 bertrand 387: (*((real8 *) (*(*(*s_etat_processus).pointeur_variable_courante)
1.1 bertrand 388: .objet).objet)) = (*point);
389: }
390: else
391: {
392: /*
393: * La variable d'intégration étant locale, l'utilisateur ne peut
394: * pas l'effacer. On ne doit donc jamais passer par ici. Si
395: * c'est néanmoins le cas, une erreur système est générée
396: * provoquant l'arrêt du programme même dans une
397: * structure IFERR.
398: */
399:
400: return;
401: }
402:
403: hauteur_pile = (*s_etat_processus).hauteur_pile_operationnelle;
404: l_position_normale = (*s_etat_processus).l_base_pile_systeme;
405:
406: (*s_etat_processus).erreur_execution = d_ex;
407: (*s_etat_processus).exception = d_ep;
408:
409: if (evaluation(s_etat_processus, s_expression, 'N') == d_erreur)
410: {
411: if ((*s_etat_processus).erreur_systeme != d_es)
412: {
413: /*
414: * Erreur système
415: */
416:
417: return;
418: }
419: else
420: {
421: /*
422: * Restauration de la pile initiale
423: */
424:
425: while((*s_etat_processus).hauteur_pile_operationnelle >
426: (unsigned long) hauteur_pile)
427: {
428: if (depilement(s_etat_processus, &((*s_etat_processus)
429: .l_base_pile), &s_objet) == d_erreur)
430: {
431: (*s_etat_processus).erreur_execution =
432: d_ex_manque_argument;
433: break;
434: }
435:
436: liberation(s_etat_processus, s_objet);
437: }
438:
439: /*
440: * Restauration de la pile système
441: */
442:
443: while((*s_etat_processus).l_base_pile_systeme !=
444: l_position_normale)
445: {
446: depilement_pile_systeme(s_etat_processus);
447:
448: if ((*s_etat_processus).erreur_systeme != d_es)
449: {
450: /*
451: * Une pile vide provoque une erreur système
452: */
453:
454: return;
455: }
456: }
457: }
458:
459: (*s_etat_processus).erreur_execution = d_ex;
460: (*s_etat_processus).exception = d_ep;
461: }
462: else
463: {
464: /*
465: * Retrait de la valeur de la pile
466: */
467:
468: if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
469: &s_objet) == d_erreur)
470: {
471: (*s_etat_processus).erreur_execution = d_ex_manque_argument;
472: return;
473: }
474:
475: if ((*s_objet).type == INT)
476: {
477: (*valeur) = (*((integer8 *) (*s_objet).objet));
478: (*validite) = d_vrai;
479: }
480: else if ((*s_objet).type == REL)
481: {
482: (*valeur) = (*((real8 *) (*s_objet).objet));
483: (*validite) = d_vrai;
484: }
485:
486: liberation(s_etat_processus, s_objet);
487: }
488:
489: return;
490: }
491:
492: // vim: ts=4