FORM v5.0.0-35-g6318119
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 MLOCK(ErrorMessageLock);
242 MesPrint("Unnormalized argument in GetFloatArgument: %a",*term,term);
243 MUNLOCK(ErrorMessageLock);
244 Terminate(-1);
245 }
246 else if ( t[0] == SYMBOL && t[1] == 4 && t[2] == PISYMBOL && t[3] == 1 ) {
247 if ( first ) {
248 mpfr_const_pi(f_out,RND);
249 first = 0;
250 }
251 else {
252 mpfr_const_pi(auxr5,RND);
253 mpfr_mul(f_out,f_out,auxr5,RND);
254 }
255 }
256 else { /* This we do not / cannot do */
257 return(-1);
258 }
259 t += t[1];
260 }
261 return(0);
262}
263
264/*
265 #] GetFloatArgument :
266 #[ GetPiArgument :
267
268 Tests for sin,cos,tan whether the argument is a simple
269 multiple of pi or can be reduced to such.
270 Return value: -1: the answer is no.
271 0-23: we have 0, 1/12*pi_,...,23/12*pi_
272 With success most values can be worked out by simple means.
273*/
274
275int GetPiArgument(PHEAD WORD *arg)
276{
277 UWORD *coef, *co, *tco, twelve, *numer, *denom, *c, *d;
278 WORD i, ii, i2, iflip, nc, nd, *t, *tn, *tstop;
279 int rem;
280/*
281 One: determine whether there is a rational coefficient and a pi_:
282*/
283 if ( *arg == -SNUMBER && arg[1] == 0 ) return(0);
284 if ( *arg == -SYMBOL && arg[1] == PISYMBOL ) return(12);
285 if ( *arg < 0 ) return(-1);
286 if ( arg[ARGHEAD] != *arg-ARGHEAD ) return(-1);
287 t = arg+ARGHEAD+1;
288 tn = arg+*arg; tstop = tn - ABS(tn[-1]);
289 if ( *t != SYMBOL || t[1] != 4 || t[2] != PISYMBOL || t[3] != 1
290 || t+t[1] != tstop ) return(-1);
291/*
292 The denominator must be a divisor of 12
293 1: copy the coefficient
294*/
295 co = coef = NumberMalloc("GetPiArgument");
296 tco = (UWORD *)tstop;
297 i = tn[-1]; if ( i < 0 ) { i = -i; iflip = 1; } else iflip = 0;
298 ii = (i-1)/2;
299 twelve = 12;
300 NCOPY(co,tco,i);
301 co = coef;
302 Mully(BHEAD co,&ii,&twelve,1);
303/*
304 Now the denominator should be 1
305*/
306 i = ii; i2 = i-1;
307 numer = co; denom = numer + i;
308 if ( i > 1 ) {
309 while ( i2 > 0 ) {
310 if ( denom[i2--] != 0 ) {
311 NumberFree(co,"GetPiArgument");
312 return(-1);
313 }
314 }
315 }
316 if ( *denom != 1 ) {
317 NumberFree(co,"GetPiArgument");
318 return(-1);
319 }
320/*
321 Now we need the numerator modulus 24.
322*/
323 if ( i == 1 ) {
324 rem = *numer % 24;
325 }
326 else {
327 c = NumberMalloc("GetPiArgument");
328 d = NumberMalloc("GetPiArgument");
329 twelve *= 2;
330 DivLong(numer,i,&twelve,1,c,&nc,d,&nd);
331 rem = *d % 24;
332 NumberFree(d,"GetPiArgument");
333 NumberFree(c,"GetPiArgument");
334 }
335 NumberFree(co,"GetPiArgument");
336 if ( iflip && rem != 0 ) rem = 24-rem;
337 return(rem);
338}
339
340/*
341 #] GetPiArgument :
342 #[ EvaluateFun :
343
344 What we need to do is:
345 1: look for a function to be treated.
346 2: make sure its argument is treatable.
347 3: call the proper mpfr_ function.
348 4: accumulate the result.
349 For some functions we have to insert 'smart' shortcuts as is
350 the case with sin_(pi_) or sin_(pi_/6) of sqrt(4/9) etc.
351 Otherwise we may have to insert a value for pi_ first.
352 There are several types of arguments:
353 a: (short) integers.
354 b: rationals.
355 c: floats.
356 We accumulate the result(s) in auxr2. The argument comes in aux5
357 and a result in auxr3 which then gets multiplied into auxr2.
358 In the end we have to combine auxr2 with whatever coefficient
359 existed already.
360
361 When the float system is started we need for aux only 3 variables
362 per thread. auxr1-auxr3. This should be done separately.
363 The main problem is the conversion of mpfr_t to float_ and/or mpf_t
364 and back.
365*/
366
367
368int EvaluateFun(PHEAD WORD *term, WORD level, WORD *pars)
369{
370 WORD *t, *tstop, *tt, *tnext, *newterm, i;
371 WORD *oldworkpointer = AT.WorkPointer, nsize, nsgn, nsgn2;
372 int retval = 0, first = 1, pimul;
373
374 tstop = term + *term; tstop -= ABS(tstop[-1]);
375 if ( AT.WorkPointer < term+*term ) AT.WorkPointer = term + *term;
376 t = term+1;
377 mpfr_set_ui(auxr2,1L,RND);
378 while ( t < tstop ) {
379 if ( pars[2] == *t ) { /* have to do this one if possible */
380TestArgument:
381/*
382 There must be a single argument, except for the AGM or atan2 functions
383*/
384 tnext = t+t[1]; tt = t+FUNHEAD; NEXTARG(tt);
385 if( *t == SYMBOL) {
386 for ( WORD* ti = t+2; ti < t+t[1]; ti+=2 ) {
387 if( ( *ti == PISYMBOL || *ti == EESYMBOL || *ti == EMSYMBOL )
388 && ( pars[2] == ALLFUNCTIONS || pars[3] == *ti ) ) {
389 if ( *ti == PISYMBOL )
390 mpfr_const_pi(auxr3,RND);
391 else if ( *ti == EESYMBOL ) {
392 mpfr_set_ui(auxr3,1,RND);
393 mpfr_exp(auxr3,auxr3,RND);
394 }
395 else if ( *ti == EMSYMBOL )
396 mpfr_const_euler(auxr3,RND);
397 if ( ti[1] != 1 )
398 mpfr_pow_si(auxr3,auxr3,ti[1],RND);
399 mpfr_mul(auxr2,auxr2,auxr3,RND);
400 ti[1] = 0;
401 first = 0;
402 }
403 }
404 goto nextfun;
405 }
406 if ( tt != tnext && *t != AGMFUNCTION && *t != ATAN2FUNCTION) goto nextfun;
407 if ( *t == SINFUNCTION ) {
408 pimul = GetPiArgument(BHEAD t+FUNHEAD);
409 if ( pimul >= 0 && pimul < 24 ) {
410 if ( pimul == 0 || pimul == 12 ) goto getout;
411 if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }
412 else nsgn = 1;
413 if ( pimul > 6 ) pimul = 12-pimul;
414 switch ( pimul ) {
415 case 1:
416label1:
417 mpfr_sqrt_ui(auxr3,3L,RND);
418 mpfr_sub_ui(auxr3,auxr3,1L,RND);
419 mpfr_sqrt_ui(auxr1,2L,RND);
420 mpfr_div_ui(auxr1,auxr1,4L,RND);
421 mpfr_mul(auxr3,auxr3,auxr1,RND);
422 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
423 mpfr_mul(auxr2,auxr2,auxr3,RND);
424 break;
425 case 2:
426label2:
427 mpfr_div_ui(auxr2,auxr2,2L,RND);
428 if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);
429 break;
430 case 3:
431label3:
432 mpfr_sqrt_ui(auxr3,2L,RND);
433 mpfr_div_ui(auxr3,auxr3,2L,RND);
434 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
435 mpfr_mul(auxr2,auxr2,auxr3,RND);
436 break;
437 case 4:
438label4:
439 mpfr_sqrt_ui(auxr3,3L,RND);
440 mpfr_div_ui(auxr3,auxr3,2L,RND);
441 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
442 mpfr_mul(auxr2,auxr2,auxr3,RND);
443 break;
444 case 5:
445label5:
446 mpfr_sqrt_ui(auxr3,3L,RND);
447 mpfr_add_ui(auxr3,auxr3,1L,RND);
448 mpfr_sqrt_ui(auxr1,2L,RND);
449 mpfr_div_ui(auxr1,auxr1,4L,RND);
450 mpfr_mul(auxr3,auxr3,auxr1,RND);
451 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
452 mpfr_mul(auxr2,auxr2,auxr3,RND);
453 break;
454 case 6:
455label6:
456 if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);
457 break;
458 }
459 *t = 0; first = 0;
460 goto nextfun;
461 }
462 }
463 else if ( *t == COSFUNCTION ) {
464 pimul = GetPiArgument(BHEAD t+FUNHEAD);
465 if ( pimul >= 0 && pimul < 24 ) {
466 if ( pimul > 12 ) pimul = 24-pimul;
467 if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -1; }
468 else nsgn = 1;
469 if ( pimul == 6 ) goto getout;
470 switch ( pimul ) {
471 case 0: goto label6;
472 case 1: goto label5;
473 case 2: goto label4;
474 case 3: goto label3;
475 case 4: goto label2;
476 case 5: goto label1;
477 }
478 }
479 }
480 else if ( *t == TANFUNCTION ) {
481 pimul = GetPiArgument(BHEAD t+FUNHEAD);
482 if ( pimul >= 0 && pimul < 24 ) {
483 if ( pimul == 6 || pimul == 18 ) goto nextfun;
484 if ( pimul == 0 || pimul == 12 ) goto getout;
485 if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }
486 else nsgn = 1;
487 if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -nsgn; }
488 switch ( pimul ) {
489 case 1:
490 mpfr_sqrt_ui(auxr3,3L,RND);
491 mpfr_sub_ui(auxr3,auxr3,2L,RND);
492 if ( nsgn > 0 ) mpfr_neg(auxr3,auxr3,RND);
493 mpfr_mul(auxr2,auxr2,auxr3,RND);
494 break;
495 case 2:
496 mpfr_sqrt_ui(auxr3,3L,RND);
497 mpfr_div_ui(auxr3,auxr3,3L,RND);
498 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
499 mpfr_mul(auxr2,auxr2,auxr3,RND);
500 break;
501 case 3:
502 break;
503 case 4:
504 mpfr_sqrt_ui(auxr3,3L,RND);
505 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
506 mpfr_mul(auxr2,auxr2,auxr3,RND);
507 break;
508 case 5:
509 mpfr_sqrt_ui(auxr3,3L,RND);
510 mpfr_add_ui(auxr3,auxr3,2L,RND);
511 if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
512 mpfr_mul(auxr2,auxr2,auxr3,RND);
513 break;
514 }
515 *t = 0; first = 0;
516 goto nextfun;
517 }
518 }
519
520 if ( *t == AGMFUNCTION || *t == ATAN2FUNCTION ) {
521 if ( GetFloatArgument(BHEAD auxr1,t,1) < 0 ) goto nextfun;
522 if ( GetFloatArgument(BHEAD auxr3,t,-2) < 0 ) goto nextfun;
523 }
524 else if ( GetFloatArgument(BHEAD auxr1,t,-1) < 0 ) goto nextfun;
525 nsgn = mpfr_sgn(auxr1);
526
527 switch ( *t ) {
528 case SQRTFUNCTION:
529 if ( nsgn < 0 ) goto nextfun;
530 if ( nsgn == 0 ) goto getout;
531 else mpfr_sqrt(auxr3,auxr1,RND);
532 mpfr_mul(auxr2,auxr2,auxr3,RND);
533 break;
534 case LNFUNCTION:
535 if ( nsgn <= 0 ) goto nextfun;
536 if ( mpfr_cmp_ui(auxr1,1L) == 0 ) goto getout;
537 else mpfr_log(auxr3,auxr1,RND);
538 mpfr_mul(auxr2,auxr2,auxr3,RND);
539 break;
540 case EXPFUNCTION:
541 mpfr_exp(auxr3,auxr1,RND);
542 mpfr_mul(auxr2,auxr2,auxr3,RND);
543 break;
544 case LI2FUNCTION: /* should be between -1 and +1 */
545 if ( nsgn == 0 ) goto getout;
546 mpfr_abs(auxr3,auxr1,RND);
547 if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
548 mpfr_li2(auxr3,auxr1,RND);
549 mpfr_mul(auxr2,auxr2,auxr3,RND);
550 break;
551 case GAMMAFUN:
552/*
553 We cannot do this when the argument is a non-positive integer
554*/
555 if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] <= 0 ) goto nextfun;
556 if ( t[FUNHEAD] == t[1]-FUNHEAD && ABS(t[t[1]-1]) ==
557 t[FUNHEAD]-ARGHEAD-1 && t[t[1]-1] < 0 ) {
558 nsize = (-t[t[1]-1]-1)/2;
559 if ( t[t[1]-1-nsize] == 1 ) {
560 for ( i = 1; i < nsize; i++ ) {
561 if ( t[t[1]-1-nsize+i] != 0 ) break;
562 }
563 if ( i >= nsize ) goto nextfun;
564 }
565 }
566 mpfr_gamma(auxr3,auxr1,RND);
567 mpfr_mul(auxr2,auxr2,auxr3,RND);
568 break;
569 case AGMFUNCTION:
570 nsgn = mpfr_sgn(auxr1);
571 nsgn2 = mpfr_sgn(auxr3);
572 if ( nsgn == 0 || nsgn2 == 0) goto getout;
573 if ( nsgn < 0 || nsgn2 < 0 ) goto nextfun;
574 mpfr_agm(auxr3,auxr1,auxr3,RND);
575 mpfr_mul(auxr2,auxr2,auxr3,RND);
576 break;
577 case SINHFUNCTION:
578 if ( nsgn == 0 ) goto getout;
579 mpfr_sinh(auxr3,auxr1,RND);
580 mpfr_mul(auxr2,auxr2,auxr3,RND);
581 break;
582 case COSHFUNCTION:
583 mpfr_cosh(auxr3,auxr1,RND);
584 mpfr_mul(auxr2,auxr2,auxr3,RND);
585 break;
586 case TANHFUNCTION:
587 if ( nsgn == 0 ) goto getout;
588 mpfr_tanh(auxr3,auxr1,RND);
589 mpfr_mul(auxr2,auxr2,auxr3,RND);
590 break;
591 case ASINHFUNCTION:
592 if ( nsgn == 0 ) goto getout;
593 mpfr_asinh(auxr3,auxr1,RND);
594 mpfr_mul(auxr2,auxr2,auxr3,RND);
595 break;
596 case ACOSHFUNCTION:
597 if ( mpfr_cmp_ui(auxr1,1L) < 0 ) goto nextfun;
598 if ( mpfr_cmp_ui(auxr1,1L) == 0 ) goto getout;
599 mpfr_acosh(auxr3,auxr1,RND);
600 mpfr_mul(auxr2,auxr2,auxr3,RND);
601 break;
602 case ATANHFUNCTION:
603 if ( nsgn == 0 ) goto getout;
604 mpfr_abs(auxr3,auxr1,RND);
605 if ( mpfr_cmp_ui(auxr3,1L) >= 0 ) goto nextfun;
606 mpfr_atanh(auxr3,auxr1,RND);
607 mpfr_mul(auxr2,auxr2,auxr3,RND);
608 break;
609 case ASINFUNCTION:
610 if ( nsgn == 0 ) goto getout;
611 mpfr_abs(auxr3,auxr1,RND);
612 if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
613 mpfr_asin(auxr3,auxr1,RND);
614 mpfr_mul(auxr2,auxr2,auxr3,RND);
615 break;
616 case ACOSFUNCTION:
617 if ( mpfr_cmp_ui(auxr1,1L) == 0 ) goto getout;
618 mpfr_abs(auxr3,auxr1,RND);
619 if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
620 mpfr_acos(auxr3,auxr1,RND);
621 mpfr_mul(auxr2,auxr2,auxr3,RND);
622 break;
623 case ATANFUNCTION:
624 if ( nsgn == 0 ) goto getout;
625 mpfr_atan(auxr3,auxr1,RND);
626 mpfr_mul(auxr2,auxr2,auxr3,RND);
627 break;
628 case ATAN2FUNCTION:
629 nsgn = mpfr_sgn(auxr1);
630 nsgn2 = mpfr_sgn(auxr3);
631 // We follow the conventions of mpfr here:
632 if ( nsgn == 0 && nsgn2 >= 0) goto getout;
633 mpfr_atan2(auxr3,auxr1,auxr3,RND);
634 mpfr_mul(auxr2,auxr2,auxr3,RND);
635 break;
636 case SINFUNCTION:
637 mpfr_sin(auxr3,auxr1,RND);
638 mpfr_mul(auxr2,auxr2,auxr3,RND);
639 break;
640 case COSFUNCTION:
641 mpfr_cos(auxr3,auxr1,RND);
642 mpfr_mul(auxr2,auxr2,auxr3,RND);
643 break;
644 case TANFUNCTION:
645 mpfr_tan(auxr3,auxr1,RND);
646 mpfr_mul(auxr2,auxr2,auxr3,RND);
647 break;
648 default:
649 goto nextfun;
650 break;
651 }
652 first = 0;
653 *t = 0;
654 goto nextfun;
655 }
656 else if ( pars[2] == ALLFUNCTIONS ) {
657/*
658 Now we have to test whether this is one of our functions
659*/
660 if ( t[1] == FUNHEAD ) goto nextfun;
661 switch ( *t ) {
662 case SQRTFUNCTION:
663 case LNFUNCTION:
664 case LI2FUNCTION:
665 case LINFUNCTION:
666 case EXPFUNCTION:
667 case ASINFUNCTION:
668 case ACOSFUNCTION:
669 case ATANFUNCTION:
670 case ATAN2FUNCTION:
671 case SINHFUNCTION:
672 case COSHFUNCTION:
673 case TANHFUNCTION:
674 case ASINHFUNCTION:
675 case ACOSHFUNCTION:
676 case ATANHFUNCTION:
677 case SINFUNCTION:
678 case COSFUNCTION:
679 case TANFUNCTION:
680 case AGMFUNCTION:
681 case SYMBOL:
682 goto TestArgument;
683 case MZV:
684 case EULER:
685 case MZVHALF:
686 default:
687 goto nextfun;
688 }
689 }
690 else goto nextfun;
691nextfun:
692 t += t[1];
693 }
694 if ( first == 1 ) return(Generator(BHEAD term,level));
695 mpfr_get_f(aux4,auxr2,RND);
696/*
697 Step 3:
698 Now the regular coefficient, if it is not 1/1.
699 We have two cases: size +- 3, or bigger.
700*/
701
702 nsize = term[*term-1];
703 if ( nsize < 0 ) { nsize = -nsize; nsgn = -1; }
704 else nsgn = 1;
705 if ( aux4->_mp_size < 0 ) {
706 aux4->_mp_size = -aux4->_mp_size;
707 nsgn = -nsgn;
708 }
709
710 if ( nsize == 3 ) {
711 if ( tstop[0] != 1 ) {
712 mpf_mul_ui(aux4,aux4,(ULONG)((UWORD)tstop[0]));
713 }
714 if ( tstop[1] != 1 ) {
715 mpf_div_ui(aux4,aux4,(ULONG)((UWORD)tstop[1]));
716 }
717 }
718 else {
719 RatToFloat(aux5,(UWORD *)tstop,nsize);
720 mpf_mul(aux4,aux4,aux5);
721 }
722/*
723 Now we have to locate possible other float_ functions.
724 Note possible incompatibilities between the mpf and mpfr formats.
725*/
726 t = term+1;
727 while ( t < tstop ) {
728 if ( *t == FLOATFUN ) {
729 UnpackFloat(aux5,t);
730 mpf_mul(aux4,aux4,aux5);
731 }
732 t += t[1];
733 }
734/*
735 Now we should compose the new term in the WorkSpace.
736*/
737 t = term+1;
738 newterm = AT.WorkPointer;
739 tt = newterm+1;
740 while ( t < tstop ) {
741 if ( *t == 0 || *t == FLOATFUN ) t += t[1];
742 else {
743 i = t[1]; NCOPY(tt,t,i);
744 }
745 }
746 PackFloat(tt,aux4);
747 tt += tt[1];
748 *tt++ = 1; *tt++ = 1; *tt++ = 3*nsgn;
749 *newterm = tt-newterm;
750 AT.WorkPointer = tt;
751 retval = Generator(BHEAD newterm,level);
752getout:
753 AT.WorkPointer = oldworkpointer;
754 return(retval);
755}
756
757/*
758 #] EvaluateFun :
759 #] mpfr_ :
760*/
Generator
int Generator(PHEAD WORD *, WORD)
Definition proces.c:3249
Generated by doxygen 1.9.8