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