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