FORM v5.0.1-23-g7a8f756
evaluate.c
Go to the documentation of this file.
1
5/* #[ License : */
6/*
7 * Copyright (C) 1984-2026 J.A.M. Vermaseren
8 * When using this file you are requested to refer to the publication
9 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
10 * This is considered a matter of courtesy as the development was paid
11 * for by FOM the Dutch physics granting agency and we would like to
12 * be able to track its scientific use to convince FOM of its value
13 * for the community.
14 *
15 * This file is part of FORM.
16 *
17 * FORM is free software: you can redistribute it and/or modify it under the
18 * terms of the GNU General Public License as published by the Free Software
19 * Foundation, either version 3 of the License, or (at your option) any later
20 * version.
21 *
22 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
23 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
24 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
25 * details.
26 *
27 * You should have received a copy of the GNU General Public License along
28 * with FORM. If not, see <http://www.gnu.org/licenses/>.
29 */
30/* #] License : */
31/*
32 #[ includes :
33*/
34
35#include "form3.h"
36#include <gmp.h>
37#include <mpfr.h>
38
39#define EXTRAPRECISION 4
40
41#define RND MPFR_RNDN
42
43#define auxr1 (((mpfr_t *)(AT.auxr_))[0])
44#define auxr2 (((mpfr_t *)(AT.auxr_))[1])
45#define auxr3 (((mpfr_t *)(AT.auxr_))[2])
46#define auxr4 (((mpfr_t *)(AT.auxr_))[3])
47#define auxr5 (((mpfr_t *)(AT.auxr_))[4])
48
49mpf_t ln2;
50
51int PackFloat(WORD *,mpf_t);
52int UnpackFloat(mpf_t,WORD *);
53void RatToFloat(mpf_t, UWORD *, int);
54void FormtoZ(mpz_t,UWORD *,WORD);
55
56/*
57 #] includes :
58 #[ mpfr_ :
59 #[ IntegerToFloatr :
60
61 Converts a Form long integer to a mpfr_ float of default size.
62 We assume that sizeof(unsigned long int) = 2*sizeof(UWORD).
63*/
64
65void IntegerToFloatr(mpfr_t result, UWORD *formlong, int longsize)
66{
67 mpz_t z;
68 mpz_init(z);
69 FormtoZ(z,formlong,longsize);
70 mpfr_set_z(result,z,RND);
71 mpz_clear(z);
72}
73
74/*
75 #] IntegerToFloatr :
76 #[ RatToFloatr :
77
78 Converts a Form rational to a gmp float of default size.
79*/
80
81void RatToFloatr(mpfr_t result, UWORD *formrat, int ratsize)
82{
83 GETIDENTITY
84 int nnum, nden;
85 UWORD *num, *den;
86 int sgn = 0;
87 if ( ratsize < 0 ) { ratsize = -ratsize; sgn = 1; }
88 nnum = nden = (ratsize-1)/2;
89 num = formrat; den = formrat+nnum;
90 while ( nnum > 1 && num[nnum-1] == 0 ) { nnum--; }
91 while ( nden > 1 && den[nden-1] == 0 ) { nden--; }
92 IntegerToFloatr(auxr4,num,nnum);
93 IntegerToFloatr(auxr5,den,nden);
94 mpfr_div(result,auxr4,auxr5,RND);
95 if ( sgn > 0 ) mpfr_neg(result,result,RND);
96}
97
98/*
99 #] RatToFloatr :
100 #[ SetfFloatPrecision :
101
102 The AC.DefaultPrecision is in bits.
103 prec is in limbs.
104 fprec is NOT in bits for mpfr_ which is slightly less than in mpf_
105 We also set the auxr1,auxr2,auxr3,auxr4,auxr5 variables.
106*/
107
108void SetfFloatPrecision(LONG prec)
109{
110/*
111 mpfr_prec_t fprec = prec * mp_bits_per_limb - EXTRAPRECISION;
112*/
113 mpfr_prec_t fprec = prec + 1;
114 mpfr_set_default_prec(fprec);
115#ifdef WITHPTHREADS
116 int totnum = AM.totalnumberofthreads, id;
117 mpfr_t *a;
118 #ifdef WITHSORTBOTS
119 totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);
120 #endif
121 for ( id = 0; id < totnum; id++ ) {
122 AB[id]->T.auxr_ = (void *)Malloc1(sizeof(mpfr_t)*5,"AB[id]->T.auxr_");
123 a = (mpfr_t *)AB[id]->T.auxr_;
124/*
125 We work here with a[0] etc because the aux1 etc contain B which
126 in the current routine would be AB[0] only
127*/
128 mpfr_inits2(fprec,a[0],a[1],a[2],a[3],a[4],(mpfr_ptr)0);
129 }
130#else
131 AT.auxr_ = (void *)Malloc1(sizeof(mpfr_t)*5,"AT.auxr_");
132 mpfr_inits2(fprec,auxr1,auxr2,auxr3,auxr4,auxr5,(mpfr_ptr)0);
133#endif
134}
135
136/*
137 #] SetfFloatPrecision :
138 #[ ClearfFloat :
139*/
140
141void ClearfFloat(void)
142{
143#ifdef WITHPTHREADS
144 int totnum = AM.totalnumberofthreads, id;
145 mpfr_t *a;
146 #ifdef WITHSORTBOTS
147 totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);
148 #endif
149 if ( AB[0]->T.auxr_ ) {
150 for ( id = 0; id < totnum; id++ ) {
151 a = (mpfr_t *)AB[id]->T.auxr_;
152 mpfr_clears(a[0],a[1],a[2],a[3],a[4],(mpfr_ptr)0);
153 M_free(AB[id]->T.auxr_,"AB[id]->T.auxr_");
154 AB[id]->T.auxr_ = 0;
155 }
156 }
157#else
158 if ( AT.auxr_ ) {
159 mpfr_clears(auxr1,auxr2,auxr3,auxr4,auxr5,(mpfr_ptr)0);
160 M_free(AT.auxr_,"AT.auxr_");
161 AT.auxr_ = 0;
162 }
163#endif
164}
165
166/*
167 #] ClearfFloat :
168 #[ GetFloatArgument :
169
170 Convert an argument to an mpfr float if possible.
171 Return value: 0: was converted successfully.
172 -1: could not convert.
173 Note: arguments with more than one term should be evaluated first
174 inside an Argument environment. If there is still more than one
175 term remaining we get the code: "could not convert".
176 par tells which argument we need. If it is negative it should also
177 be the last argument.
178*/
179
180int GetFloatArgument(PHEAD mpfr_t f_out,WORD *fun,int par)
181{
182 WORD *term, *tn, *t, *tstop, first, ncoef, *arg, *argp, abspar;
183 arg = fun+FUNHEAD;
184 abspar = ABS(par);
185 while ( arg < fun+fun[1] && abspar > 1 ) { abspar--; NEXTARG(arg) }
186 if ( arg >= fun+fun[1] || abspar!= 1 ) return(-1);
187 if ( par < 0 ) {
188 argp = arg; NEXTARG(argp); if ( argp != fun+fun[1] ) return(-1);
189 }
190 if ( *arg < 0 ) {
191 if ( *arg == -SNUMBER ) {
192 mpfr_set_si(f_out,(LONG)(arg[1]),RND);
193 }
194 else if ( *arg == -SYMBOL && ( arg[1] == PISYMBOL ) ) {
195 mpfr_const_pi(f_out,RND);
196 }
197 else if ( *arg == -INDEX && arg[1] >= 0 && arg[1] < AM.OffsetIndex ) {
198 mpfr_set_ui(f_out,(ULONG)(arg[1]),RND);
199 }
200 else { return(-1); }
201 return(0);
202 }
203 else if ( arg[0] != arg[ARGHEAD]+ARGHEAD ) { /* more than one term */
204 return(-1);
205 }
206 term = arg+ARGHEAD;
207 tn = term+*term;
208 tstop = tn-ABS(tn[-1]);
209 t = term+1;
210 first = 1;
211 while ( t <= tstop ) {
212 if ( t == tstop ) { /* Fraction */
213 if ( first ) RatToFloatr(f_out,(UWORD *)tstop,tn[-1]);
214 else {
215 RatToFloatr(auxr4,(UWORD *)tstop,tn[-1]);
216 mpfr_mul(f_out,f_out,auxr4,RND);
217 }
218 return(0);
219 }
220 else if ( *t == FLOATFUN ) {
221 if ( t+t[1] != tstop ) return(-1);
222 if ( UnpackFloat(aux5,t) < 0 ) return(-1);
223 if ( first ) {
224 mpfr_set_f(f_out,aux5,RND);
225 first = 0;
226 }
227 else {
228 mpfr_set_f(auxr4,aux5,RND);
229 mpfr_mul(f_out,f_out,auxr4,RND);
230 }
231 first = 0;
232 if ( tn[-1] < 0 ) { /* change sign */
233 ncoef = -tn[-1];
234 mpfr_neg(f_out,f_out,RND);
235 }
236 else ncoef = tn[-1];
237 if ( ncoef == 3 && tn[-2] == 1 && tn[-3] == 1 ) return(0);
238/*
239 If the argument was properly normalized we are not supposed to come here.
240*/
241/* INTERNAL_ERROR_EXCL_START */
242 MLOCK(ErrorMessageLock);
243 MesPrint("!>Unnormalized argument in GetFloatArgument: %a",*term,term);
244 MUNLOCK(ErrorMessageLock);
245 Terminate(-1);
246/* INTERNAL_ERROR_EXCL_STOP */
247 }
248 else if ( t[0] == SYMBOL && t[1] == 4 && t[2] == PISYMBOL && t[3] == 1 ) {
249 if ( first ) {
250 mpfr_const_pi(f_out,RND);
251 first = 0;
252 }
253 else {
254 mpfr_const_pi(auxr5,RND);
255 mpfr_mul(f_out,f_out,auxr5,RND);
256 }
257 }
258 else { /* This we do not / cannot do */
259 return(-1);
260 }
261 t += t[1];
262 }
263 return(0);
264}
265
266/*
267 #] GetFloatArgument :
268 #[ GetPiArgument :
269
270 Tests for sin,cos,tan whether the argument is a simple
271 multiple of pi or can be reduced to such.
272 Return value: -1: the answer is no.
273 0-23: we have 0, 1/12*pi_,...,23/12*pi_
274 With success most values can be worked out by simple means.
275*/
276
277int GetPiArgument(PHEAD WORD *arg)
278{
279 UWORD *coef, *co, *tco, twelve, *numer, *denom, *c, *d;
280 WORD i, ii, i2, iflip, nc, nd, *t, *tn, *tstop;
281 int rem;
282/*
283 One: determine whether there is a rational coefficient and a pi_:
284*/
285 if ( *arg == -SNUMBER && arg[1] == 0 ) return(0);
286 if ( *arg == -SYMBOL && arg[1] == PISYMBOL ) return(12);
287 if ( *arg < 0 ) return(-1);
288 if ( arg[ARGHEAD] != *arg-ARGHEAD ) return(-1);
289 t = arg+ARGHEAD+1;
290 tn = arg+*arg; tstop = tn - ABS(tn[-1]);
291 if ( *t != SYMBOL || t[1] != 4 || t[2] != PISYMBOL || t[3] != 1
292 || t+t[1] != tstop ) return(-1);
293/*
294 The denominator must be a divisor of 12
295 1: copy the coefficient
296*/
297 co = coef = NumberMalloc("GetPiArgument");
298 tco = (UWORD *)tstop;
299 i = tn[-1]; if ( i < 0 ) { i = -i; iflip = 1; } else iflip = 0;
300 ii = (i-1)/2;
301 twelve = 12;
302 NCOPY(co,tco,i);
303 co = coef;
304 Mully(BHEAD co,&ii,&twelve,1);
305/*
306 Now the denominator should be 1
307*/
308 i = ii; i2 = i-1;
309 numer = co; denom = numer + i;
310 if ( i > 1 ) {
311 while ( i2 > 0 ) {
312 if ( denom[i2--] != 0 ) {
313 NumberFree(co,"GetPiArgument");
314 return(-1);
315 }
316 }
317 }
318 if ( *denom != 1 ) {
319 NumberFree(co,"GetPiArgument");
320 return(-1);
321 }
322/*
323 Now we need the numerator modulus 24.
324*/
325 if ( i == 1 ) {
326 rem = *numer % 24;
327 }
328 else {
329 c = NumberMalloc("GetPiArgument");
330 d = NumberMalloc("GetPiArgument");
331 twelve *= 2;
332 DivLong(numer,i,&twelve,1,c,&nc,d,&nd);
333 rem = *d % 24;
334 NumberFree(d,"GetPiArgument");
335 NumberFree(c,"GetPiArgument");
336 }
337 NumberFree(co,"GetPiArgument");
338 if ( iflip && rem != 0 ) rem = 24-rem;
339 return(rem);
340}
341
342/*
343 #] GetPiArgument :
344 #[ EvaluateFun :
345
346 What we need to do is:
347 1: look for a function to be treated.
348 2: make sure its argument is treatable.
349 3: call the proper mpfr_ function.
350 4: accumulate the result.
351 For some functions we have to insert 'smart' shortcuts as is
352 the case with sin_(pi_) or sin_(pi_/6) of sqrt(4/9) etc.
353 Otherwise we may have to insert a value for pi_ first.
354 There are several types of arguments:
355 a: (short) integers.
356 b: rationals.
357 c: floats.
358 We accumulate the result(s) in auxr2. The argument comes in aux5
359 and a result in auxr3 which then gets multiplied into auxr2.
360 In the end we have to combine auxr2 with whatever coefficient
361 existed already.
362
363 When the float system is started we need for aux only 3 variables
364 per thread. auxr1-auxr3. This should be done separately.
365 The main problem is the conversion of mpfr_t to float_ and/or mpf_t
366 and back.
367*/
368
369
370int EvaluateFun(PHEAD WORD *term, WORD level, WORD *pars)
371{
372 WORD *t, *tstop, *tt, *tnext, *newterm, i;
373 WORD *oldworkpointer = AT.WorkPointer, nsize, nsgn, nsgn2;
374 int retval = 0, first = 1, pimul;
375
376 tstop = term + *term; tstop -= ABS(tstop[-1]);
377 if ( AT.WorkPointer < term+*term ) AT.WorkPointer = term + *term;
378 t = term+1;
379 mpfr_set_ui(auxr2,1L,RND);
380 while ( t < tstop ) {
381 if ( pars[2] == *t ) { /* have to do this one if possible */
382TestArgument:
383/*
384 There must be a single argument, except for the AGM or atan2 functions
385*/
386 tnext = t+t[1]; tt = t+FUNHEAD; NEXTARG(tt);
387 if( *t == SYMBOL) {
388 for ( WORD* ti = t+2; ti < t+t[1]; ti+=2 ) {
389 if( ( *ti == PISYMBOL || *ti == EESYMBOL || *ti == EMSYMBOL )
390 && ( pars[2] == ALLFUNCTIONS || pars[3] == *ti ) ) {
391 if ( *ti == PISYMBOL )
392 mpfr_const_pi(auxr3,RND);
393 else if ( *ti == EESYMBOL ) {
394 mpfr_set_ui(auxr3,1,RND);
395 mpfr_exp(auxr3,auxr3,RND);
396 }
397 else if ( *ti == EMSYMBOL )
398 mpfr_const_euler(auxr3,RND);
399 if ( ti[1] != 1 )
400 mpfr_pow_si(auxr3,auxr3,ti[1],RND);
401 mpfr_mul(auxr2,auxr2,auxr3,RND);
402 ti[1] = 0;
403 first = 0;
404 }
405 }
406 goto nextfun;
407 }
408 if ( tt != tnext && *t != AGMFUNCTION && *t != ATAN2FUNCTION) goto nextfun;
409 if ( *t == SINFUNCTION ) {
410 pimul = GetPiArgument(BHEAD t+FUNHEAD);
411 if ( pimul >= 0 && pimul < 24 ) {
412 if ( pimul == 0 || pimul == 12 ) goto getout;
413 if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }
414 else nsgn = 1;
415 if ( pimul > 6 ) pimul = 12-pimul;
416 switch ( pimul ) {
417 case 1:
418label1:
419 mpfr_sqrt_ui(auxr3,3L,RND);
420 mpfr_sub_ui(auxr3,auxr3,1L,RND);
421 mpfr_sqrt_ui(auxr1,2L,RND);
422 mpfr_div_ui(auxr1,auxr1,4L,RND);
423 mpfr_mul(auxr3,auxr3,auxr1,RND);
424 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
425 mpfr_mul(auxr2,auxr2,auxr3,RND);
426 break;
427 case 2:
428label2:
429 mpfr_div_ui(auxr2,auxr2,2L,RND);
430 if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);
431 break;
432 case 3:
433label3:
434 mpfr_sqrt_ui(auxr3,2L,RND);
435 mpfr_div_ui(auxr3,auxr3,2L,RND);
436 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
437 mpfr_mul(auxr2,auxr2,auxr3,RND);
438 break;
439 case 4:
440label4:
441 mpfr_sqrt_ui(auxr3,3L,RND);
442 mpfr_div_ui(auxr3,auxr3,2L,RND);
443 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
444 mpfr_mul(auxr2,auxr2,auxr3,RND);
445 break;
446 case 5:
447label5:
448 mpfr_sqrt_ui(auxr3,3L,RND);
449 mpfr_add_ui(auxr3,auxr3,1L,RND);
450 mpfr_sqrt_ui(auxr1,2L,RND);
451 mpfr_div_ui(auxr1,auxr1,4L,RND);
452 mpfr_mul(auxr3,auxr3,auxr1,RND);
453 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
454 mpfr_mul(auxr2,auxr2,auxr3,RND);
455 break;
456 case 6:
457label6:
458 if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);
459 break;
460 }
461 *t = 0; first = 0;
462 goto nextfun;
463 }
464 }
465 else if ( *t == COSFUNCTION ) {
466 pimul = GetPiArgument(BHEAD t+FUNHEAD);
467 if ( pimul >= 0 && pimul < 24 ) {
468 if ( pimul > 12 ) pimul = 24-pimul;
469 if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -1; }
470 else nsgn = 1;
471 if ( pimul == 6 ) goto getout;
472 switch ( pimul ) {
473 case 0: goto label6;
474 case 1: goto label5;
475 case 2: goto label4;
476 case 3: goto label3;
477 case 4: goto label2;
478 case 5: goto label1;
479 }
480 }
481 }
482 else if ( *t == TANFUNCTION ) {
483 pimul = GetPiArgument(BHEAD t+FUNHEAD);
484 if ( pimul >= 0 && pimul < 24 ) {
485 if ( pimul == 6 || pimul == 18 ) goto nextfun;
486 if ( pimul == 0 || pimul == 12 ) goto getout;
487 if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }
488 else nsgn = 1;
489 if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -nsgn; }
490 switch ( pimul ) {
491 case 1:
492 mpfr_sqrt_ui(auxr3,3L,RND);
493 mpfr_sub_ui(auxr3,auxr3,2L,RND);
494 if ( nsgn > 0 ) mpfr_neg(auxr3,auxr3,RND);
495 mpfr_mul(auxr2,auxr2,auxr3,RND);
496 break;
497 case 2:
498 mpfr_sqrt_ui(auxr3,3L,RND);
499 mpfr_div_ui(auxr3,auxr3,3L,RND);
500 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
501 mpfr_mul(auxr2,auxr2,auxr3,RND);
502 break;
503 case 3:
504 break;
505 case 4:
506 mpfr_sqrt_ui(auxr3,3L,RND);
507 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
508 mpfr_mul(auxr2,auxr2,auxr3,RND);
509 break;
510 case 5:
511 mpfr_sqrt_ui(auxr3,3L,RND);
512 mpfr_add_ui(auxr3,auxr3,2L,RND);
513 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
514 mpfr_mul(auxr2,auxr2,auxr3,RND);
515 break;
516 }
517 *t = 0; first = 0;
518 goto nextfun;
519 }
520 }
521
522 if ( *t == AGMFUNCTION || *t == ATAN2FUNCTION ) {
523 if ( GetFloatArgument(BHEAD auxr1,t,1) < 0 ) goto nextfun;
524 if ( GetFloatArgument(BHEAD auxr3,t,-2) < 0 ) goto nextfun;
525 }
526 else if ( GetFloatArgument(BHEAD auxr1,t,-1) < 0 ) goto nextfun;
527 nsgn = mpfr_sgn(auxr1);
528
529 switch ( *t ) {
530 case SQRTFUNCTION:
531 if ( nsgn < 0 ) goto nextfun;
532 if ( nsgn == 0 ) goto getout;
533 else mpfr_sqrt(auxr3,auxr1,RND);
534 mpfr_mul(auxr2,auxr2,auxr3,RND);
535 break;
536 case LNFUNCTION:
537 if ( nsgn <= 0 ) goto nextfun;
538 if ( mpfr_cmp_ui(auxr1,1L) == 0 ) goto getout;
539 else mpfr_log(auxr3,auxr1,RND);
540 mpfr_mul(auxr2,auxr2,auxr3,RND);
541 break;
542 case EXPFUNCTION:
543 mpfr_exp(auxr3,auxr1,RND);
544 mpfr_mul(auxr2,auxr2,auxr3,RND);
545 break;
546 case LI2FUNCTION: /* should be between -1 and +1 */
547 if ( nsgn == 0 ) goto getout;
548 mpfr_abs(auxr3,auxr1,RND);
549 if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
550 mpfr_li2(auxr3,auxr1,RND);
551 mpfr_mul(auxr2,auxr2,auxr3,RND);
552 break;
553 case GAMMAFUN:
554/*
555 We cannot do this when the argument is a non-positive integer
556*/
557 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] <= 0 ) goto nextfun;
558 if ( t[FUNHEAD] == t[1]-FUNHEAD && ABS(t[t[1]-1]) ==
559 t[FUNHEAD]-ARGHEAD-1 && t[t[1]-1] < 0 ) {
560 nsize = (-t[t[1]-1]-1)/2;
561 if ( t[t[1]-1-nsize] == 1 ) {
562 for ( i = 1; i < nsize; i++ ) {
563 if ( t[t[1]-1-nsize+i] != 0 ) break;
564 }
565 if ( i >= nsize ) goto nextfun;
566 }
567 }
568 mpfr_gamma(auxr3,auxr1,RND);
569 mpfr_mul(auxr2,auxr2,auxr3,RND);
570 break;
571 case AGMFUNCTION:
572 nsgn = mpfr_sgn(auxr1);
573 nsgn2 = mpfr_sgn(auxr3);
574 if ( nsgn == 0 || nsgn2 == 0) goto getout;
575 if ( nsgn < 0 || nsgn2 < 0 ) goto nextfun;
576 mpfr_agm(auxr3,auxr1,auxr3,RND);
577 mpfr_mul(auxr2,auxr2,auxr3,RND);
578 break;
579 case SINHFUNCTION:
580 if ( nsgn == 0 ) goto getout;
581 mpfr_sinh(auxr3,auxr1,RND);
582 mpfr_mul(auxr2,auxr2,auxr3,RND);
583 break;
584 case COSHFUNCTION:
585 mpfr_cosh(auxr3,auxr1,RND);
586 mpfr_mul(auxr2,auxr2,auxr3,RND);
587 break;
588 case TANHFUNCTION:
589 if ( nsgn == 0 ) goto getout;
590 mpfr_tanh(auxr3,auxr1,RND);
591 mpfr_mul(auxr2,auxr2,auxr3,RND);
592 break;
593 case ASINHFUNCTION:
594 if ( nsgn == 0 ) goto getout;
595 mpfr_asinh(auxr3,auxr1,RND);
596 mpfr_mul(auxr2,auxr2,auxr3,RND);
597 break;
598 case ACOSHFUNCTION:
599 if ( mpfr_cmp_ui(auxr1,1L) < 0 ) goto nextfun;
600 if ( mpfr_cmp_ui(auxr1,1L) == 0 ) goto getout;
601 mpfr_acosh(auxr3,auxr1,RND);
602 mpfr_mul(auxr2,auxr2,auxr3,RND);
603 break;
604 case ATANHFUNCTION:
605 if ( nsgn == 0 ) goto getout;
606 mpfr_abs(auxr3,auxr1,RND);
607 if ( mpfr_cmp_ui(auxr3,1L) >= 0 ) goto nextfun;
608 mpfr_atanh(auxr3,auxr1,RND);
609 mpfr_mul(auxr2,auxr2,auxr3,RND);
610 break;
611 case ASINFUNCTION:
612 if ( nsgn == 0 ) goto getout;
613 mpfr_abs(auxr3,auxr1,RND);
614 if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
615 mpfr_asin(auxr3,auxr1,RND);
616 mpfr_mul(auxr2,auxr2,auxr3,RND);
617 break;
618 case ACOSFUNCTION:
619 if ( mpfr_cmp_ui(auxr1,1L) == 0 ) goto getout;
620 mpfr_abs(auxr3,auxr1,RND);
621 if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
622 mpfr_acos(auxr3,auxr1,RND);
623 mpfr_mul(auxr2,auxr2,auxr3,RND);
624 break;
625 case ATANFUNCTION:
626 if ( nsgn == 0 ) goto getout;
627 mpfr_atan(auxr3,auxr1,RND);
628 mpfr_mul(auxr2,auxr2,auxr3,RND);
629 break;
630 case ATAN2FUNCTION:
631 nsgn = mpfr_sgn(auxr1);
632 nsgn2 = mpfr_sgn(auxr3);
633 // We follow the conventions of mpfr here:
634 if ( nsgn == 0 && nsgn2 >= 0) goto getout;
635 mpfr_atan2(auxr3,auxr1,auxr3,RND);
636 mpfr_mul(auxr2,auxr2,auxr3,RND);
637 break;
638 case SINFUNCTION:
639 mpfr_sin(auxr3,auxr1,RND);
640 mpfr_mul(auxr2,auxr2,auxr3,RND);
641 break;
642 case COSFUNCTION:
643 mpfr_cos(auxr3,auxr1,RND);
644 mpfr_mul(auxr2,auxr2,auxr3,RND);
645 break;
646 case TANFUNCTION:
647 mpfr_tan(auxr3,auxr1,RND);
648 mpfr_mul(auxr2,auxr2,auxr3,RND);
649 break;
650 default:
651 goto nextfun;
652 break;
653 }
654 first = 0;
655 *t = 0;
656 goto nextfun;
657 }
658 else if ( pars[2] == ALLFUNCTIONS ) {
659/*
660 Now we have to test whether this is one of our functions
661*/
662 if ( t[1] == FUNHEAD ) goto nextfun;
663 switch ( *t ) {
664 case SQRTFUNCTION:
665 case LNFUNCTION:
666 case LI2FUNCTION:
667 case LINFUNCTION:
668 case EXPFUNCTION:
669 case ASINFUNCTION:
670 case ACOSFUNCTION:
671 case ATANFUNCTION:
672 case ATAN2FUNCTION:
673 case SINHFUNCTION:
674 case COSHFUNCTION:
675 case TANHFUNCTION:
676 case ASINHFUNCTION:
677 case ACOSHFUNCTION:
678 case ATANHFUNCTION:
679 case SINFUNCTION:
680 case COSFUNCTION:
681 case TANFUNCTION:
682 case AGMFUNCTION:
683 case SYMBOL:
684 goto TestArgument;
685 case MZV:
686 case EULER:
687 case MZVHALF:
688 default:
689 goto nextfun;
690 }
691 }
692 else goto nextfun;
693nextfun:
694 t += t[1];
695 }
696 if ( first == 1 ) return(Generator(BHEAD term,level));
697 mpfr_get_f(aux4,auxr2,RND);
698/*
699 Step 3:
700 Now the regular coefficient, if it is not 1/1.
701 We have two cases: size +- 3, or bigger.
702*/
703
704 nsize = term[*term-1];
705 if ( nsize < 0 ) { nsize = -nsize; nsgn = -1; }
706 else nsgn = 1;
707 if ( aux4->_mp_size < 0 ) {
708 aux4->_mp_size = -aux4->_mp_size;
709 nsgn = -nsgn;
710 }
711
712 if ( nsize == 3 ) {
713 if ( tstop[0] != 1 ) {
714 mpf_mul_ui(aux4,aux4,(ULONG)((UWORD)tstop[0]));
715 }
716 if ( tstop[1] != 1 ) {
717 mpf_div_ui(aux4,aux4,(ULONG)((UWORD)tstop[1]));
718 }
719 }
720 else {
721 RatToFloat(aux5,(UWORD *)tstop,nsize);
722 mpf_mul(aux4,aux4,aux5);
723 }
724/*
725 Now we have to locate possible other float_ functions.
726 Note possible incompatibilities between the mpf and mpfr formats.
727*/
728 t = term+1;
729 while ( t < tstop ) {
730 if ( *t == FLOATFUN ) {
731 UnpackFloat(aux5,t);
732 mpf_mul(aux4,aux4,aux5);
733 }
734 t += t[1];
735 }
736/*
737 Now we should compose the new term in the WorkSpace.
738*/
739 t = term+1;
740 newterm = AT.WorkPointer;
741 tt = newterm+1;
742 while ( t < tstop ) {
743 if ( *t == 0 || *t == FLOATFUN ) t += t[1];
744 else {
745 i = t[1]; NCOPY(tt,t,i);
746 }
747 }
748 PackFloat(tt,aux4);
749 tt += tt[1];
750 *tt++ = 1; *tt++ = 1; *tt++ = 3*nsgn;
751 *newterm = tt-newterm;
752 AT.WorkPointer = tt;
753 retval = Generator(BHEAD newterm,level);
754getout:
755 AT.WorkPointer = oldworkpointer;
756 return(retval);
757}
758
759/*
760 #] EvaluateFun :
761 #] mpfr_ :
762*/
Generator
int Generator(PHEAD WORD *, WORD)
Definition proces.c:3255
Generated by doxygen 1.9.8