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