OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgfortran/] [intrinsics/] [c99_functions.c] - Blame information for rev 733

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 733 jeremybenn
/* Implementation of various C99 functions
2
   Copyright (C) 2004, 2009, 2010 Free Software Foundation, Inc.
3
 
4
This file is part of the GNU Fortran 95 runtime library (libgfortran).
5
 
6
Libgfortran is free software; you can redistribute it and/or
7
modify it under the terms of the GNU General Public
8
License as published by the Free Software Foundation; either
9
version 3 of the License, or (at your option) any later version.
10
 
11
Libgfortran is distributed in the hope that it will be useful,
12
but WITHOUT ANY WARRANTY; without even the implied warranty of
13
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14
GNU General Public License for more details.
15
 
16
Under Section 7 of GPL version 3, you are granted additional
17
permissions described in the GCC Runtime Library Exception, version
18
3.1, as published by the Free Software Foundation.
19
 
20
You should have received a copy of the GNU General Public License and
21
a copy of the GCC Runtime Library Exception along with this program;
22
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23
<http://www.gnu.org/licenses/>.  */
24
 
25
#include "config.h"
26
 
27
#define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW
28
#include "libgfortran.h"
29
 
30
/* IRIX's <math.h> declares a non-C99 compliant implementation of cabs,
31
   which takes two floating point arguments instead of a single complex.
32
   If <complex.h> is missing this prevents building of c99_functions.c.
33
   To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}.  */
34
 
35
#if defined(__sgi__) && !defined(HAVE_COMPLEX_H)
36
#undef HAVE_CABS
37
#undef HAVE_CABSF
38
#undef HAVE_CABSL
39
#define cabs __gfc_cabs
40
#define cabsf __gfc_cabsf
41
#define cabsl __gfc_cabsl
42
#endif
43
 
44
/* Tru64's <math.h> declares a non-C99 compliant implementation of cabs,
45
   which takes two floating point arguments instead of a single complex.
46
   To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}.  */
47
 
48
#ifdef __osf__
49
#undef HAVE_CABS
50
#undef HAVE_CABSF
51
#undef HAVE_CABSL
52
#define cabs __gfc_cabs
53
#define cabsf __gfc_cabsf
54
#define cabsl __gfc_cabsl
55
#endif
56
 
57
/* On a C99 system "I" (with I*I = -1) should be defined in complex.h;
58
   if not, we define a fallback version here.  */
59
#ifndef I
60
# if defined(_Imaginary_I)
61
#   define I _Imaginary_I
62
# elif defined(_Complex_I)
63
#   define I _Complex_I
64
# else
65
#   define I (1.0fi)
66
# endif
67
#endif
68
 
69
/* Prototypes are included to silence -Wstrict-prototypes
70
   -Wmissing-prototypes.  */
71
 
72
 
73
/* Wrappers for systems without the various C99 single precision Bessel
74
   functions.  */
75
 
76
#if defined(HAVE_J0) && ! defined(HAVE_J0F)
77
#define HAVE_J0F 1
78
float j0f (float);
79
 
80
float
81
j0f (float x)
82
{
83
  return (float) j0 ((double) x);
84
}
85
#endif
86
 
87
#if defined(HAVE_J1) && !defined(HAVE_J1F)
88
#define HAVE_J1F 1
89
float j1f (float);
90
 
91
float j1f (float x)
92
{
93
  return (float) j1 ((double) x);
94
}
95
#endif
96
 
97
#if defined(HAVE_JN) && !defined(HAVE_JNF)
98
#define HAVE_JNF 1
99
float jnf (int, float);
100
 
101
float
102
jnf (int n, float x)
103
{
104
  return (float) jn (n, (double) x);
105
}
106
#endif
107
 
108
#if defined(HAVE_Y0) && !defined(HAVE_Y0F)
109
#define HAVE_Y0F 1
110
float y0f (float);
111
 
112
float
113
y0f (float x)
114
{
115
  return (float) y0 ((double) x);
116
}
117
#endif
118
 
119
#if defined(HAVE_Y1) && !defined(HAVE_Y1F)
120
#define HAVE_Y1F 1
121
float y1f (float);
122
 
123
float
124
y1f (float x)
125
{
126
  return (float) y1 ((double) x);
127
}
128
#endif
129
 
130
#if defined(HAVE_YN) && !defined(HAVE_YNF)
131
#define HAVE_YNF 1
132
float ynf (int, float);
133
 
134
float
135
ynf (int n, float x)
136
{
137
  return (float) yn (n, (double) x);
138
}
139
#endif
140
 
141
 
142
/* Wrappers for systems without the C99 erff() and erfcf() functions.  */
143
 
144
#if defined(HAVE_ERF) && !defined(HAVE_ERFF)
145
#define HAVE_ERFF 1
146
float erff (float);
147
 
148
float
149
erff (float x)
150
{
151
  return (float) erf ((double) x);
152
}
153
#endif
154
 
155
#if defined(HAVE_ERFC) && !defined(HAVE_ERFCF)
156
#define HAVE_ERFCF 1
157
float erfcf (float);
158
 
159
float
160
erfcf (float x)
161
{
162
  return (float) erfc ((double) x);
163
}
164
#endif
165
 
166
 
167
#ifndef HAVE_ACOSF
168
#define HAVE_ACOSF 1
169
float acosf (float x);
170
 
171
float
172
acosf (float x)
173
{
174
  return (float) acos (x);
175
}
176
#endif
177
 
178
#if HAVE_ACOSH && !HAVE_ACOSHF
179
float acoshf (float x);
180
 
181
float
182
acoshf (float x)
183
{
184
  return (float) acosh ((double) x);
185
}
186
#endif
187
 
188
#ifndef HAVE_ASINF
189
#define HAVE_ASINF 1
190
float asinf (float x);
191
 
192
float
193
asinf (float x)
194
{
195
  return (float) asin (x);
196
}
197
#endif
198
 
199
#if HAVE_ASINH && !HAVE_ASINHF
200
float asinhf (float x);
201
 
202
float
203
asinhf (float x)
204
{
205
  return (float) asinh ((double) x);
206
}
207
#endif
208
 
209
#ifndef HAVE_ATAN2F
210
#define HAVE_ATAN2F 1
211
float atan2f (float y, float x);
212
 
213
float
214
atan2f (float y, float x)
215
{
216
  return (float) atan2 (y, x);
217
}
218
#endif
219
 
220
#ifndef HAVE_ATANF
221
#define HAVE_ATANF 1
222
float atanf (float x);
223
 
224
float
225
atanf (float x)
226
{
227
  return (float) atan (x);
228
}
229
#endif
230
 
231
#if HAVE_ATANH && !HAVE_ATANHF
232
float atanhf (float x);
233
 
234
float
235
atanhf (float x)
236
{
237
  return (float) atanh ((double) x);
238
}
239
#endif
240
 
241
#ifndef HAVE_CEILF
242
#define HAVE_CEILF 1
243
float ceilf (float x);
244
 
245
float
246
ceilf (float x)
247
{
248
  return (float) ceil (x);
249
}
250
#endif
251
 
252
#ifndef HAVE_COPYSIGNF
253
#define HAVE_COPYSIGNF 1
254
float copysignf (float x, float y);
255
 
256
float
257
copysignf (float x, float y)
258
{
259
  return (float) copysign (x, y);
260
}
261
#endif
262
 
263
#ifndef HAVE_COSF
264
#define HAVE_COSF 1
265
float cosf (float x);
266
 
267
float
268
cosf (float x)
269
{
270
  return (float) cos (x);
271
}
272
#endif
273
 
274
#ifndef HAVE_COSHF
275
#define HAVE_COSHF 1
276
float coshf (float x);
277
 
278
float
279
coshf (float x)
280
{
281
  return (float) cosh (x);
282
}
283
#endif
284
 
285
#ifndef HAVE_EXPF
286
#define HAVE_EXPF 1
287
float expf (float x);
288
 
289
float
290
expf (float x)
291
{
292
  return (float) exp (x);
293
}
294
#endif
295
 
296
#ifndef HAVE_FABSF
297
#define HAVE_FABSF 1
298
float fabsf (float x);
299
 
300
float
301
fabsf (float x)
302
{
303
  return (float) fabs (x);
304
}
305
#endif
306
 
307
#ifndef HAVE_FLOORF
308
#define HAVE_FLOORF 1
309
float floorf (float x);
310
 
311
float
312
floorf (float x)
313
{
314
  return (float) floor (x);
315
}
316
#endif
317
 
318
#ifndef HAVE_FMODF
319
#define HAVE_FMODF 1
320
float fmodf (float x, float y);
321
 
322
float
323
fmodf (float x, float y)
324
{
325
  return (float) fmod (x, y);
326
}
327
#endif
328
 
329
#ifndef HAVE_FREXPF
330
#define HAVE_FREXPF 1
331
float frexpf (float x, int *exp);
332
 
333
float
334
frexpf (float x, int *exp)
335
{
336
  return (float) frexp (x, exp);
337
}
338
#endif
339
 
340
#ifndef HAVE_HYPOTF
341
#define HAVE_HYPOTF 1
342
float hypotf (float x, float y);
343
 
344
float
345
hypotf (float x, float y)
346
{
347
  return (float) hypot (x, y);
348
}
349
#endif
350
 
351
#ifndef HAVE_LOGF
352
#define HAVE_LOGF 1
353
float logf (float x);
354
 
355
float
356
logf (float x)
357
{
358
  return (float) log (x);
359
}
360
#endif
361
 
362
#ifndef HAVE_LOG10F
363
#define HAVE_LOG10F 1
364
float log10f (float x);
365
 
366
float
367
log10f (float x)
368
{
369
  return (float) log10 (x);
370
}
371
#endif
372
 
373
#ifndef HAVE_SCALBN
374
#define HAVE_SCALBN 1
375
double scalbn (double x, int y);
376
 
377
double
378
scalbn (double x, int y)
379
{
380
#if (FLT_RADIX == 2) && defined(HAVE_LDEXP)
381
  return ldexp (x, y);
382
#else
383
  return x * pow (FLT_RADIX, y);
384
#endif
385
}
386
#endif
387
 
388
#ifndef HAVE_SCALBNF
389
#define HAVE_SCALBNF 1
390
float scalbnf (float x, int y);
391
 
392
float
393
scalbnf (float x, int y)
394
{
395
  return (float) scalbn (x, y);
396
}
397
#endif
398
 
399
#ifndef HAVE_SINF
400
#define HAVE_SINF 1
401
float sinf (float x);
402
 
403
float
404
sinf (float x)
405
{
406
  return (float) sin (x);
407
}
408
#endif
409
 
410
#ifndef HAVE_SINHF
411
#define HAVE_SINHF 1
412
float sinhf (float x);
413
 
414
float
415
sinhf (float x)
416
{
417
  return (float) sinh (x);
418
}
419
#endif
420
 
421
#ifndef HAVE_SQRTF
422
#define HAVE_SQRTF 1
423
float sqrtf (float x);
424
 
425
float
426
sqrtf (float x)
427
{
428
  return (float) sqrt (x);
429
}
430
#endif
431
 
432
#ifndef HAVE_TANF
433
#define HAVE_TANF 1
434
float tanf (float x);
435
 
436
float
437
tanf (float x)
438
{
439
  return (float) tan (x);
440
}
441
#endif
442
 
443
#ifndef HAVE_TANHF
444
#define HAVE_TANHF 1
445
float tanhf (float x);
446
 
447
float
448
tanhf (float x)
449
{
450
  return (float) tanh (x);
451
}
452
#endif
453
 
454
#ifndef HAVE_TRUNC
455
#define HAVE_TRUNC 1
456
double trunc (double x);
457
 
458
double
459
trunc (double x)
460
{
461
  if (!isfinite (x))
462
    return x;
463
 
464
  if (x < 0.0)
465
    return - floor (-x);
466
  else
467
    return floor (x);
468
}
469
#endif
470
 
471
#ifndef HAVE_TRUNCF
472
#define HAVE_TRUNCF 1
473
float truncf (float x);
474
 
475
float
476
truncf (float x)
477
{
478
  return (float) trunc (x);
479
}
480
#endif
481
 
482
#ifndef HAVE_NEXTAFTERF
483
#define HAVE_NEXTAFTERF 1
484
/* This is a portable implementation of nextafterf that is intended to be
485
   independent of the floating point format or its in memory representation.
486
   This implementation works correctly with denormalized values.  */
487
float nextafterf (float x, float y);
488
 
489
float
490
nextafterf (float x, float y)
491
{
492
  /* This variable is marked volatile to avoid excess precision problems
493
     on some platforms, including IA-32.  */
494
  volatile float delta;
495
  float absx, denorm_min;
496
 
497
  if (isnan (x) || isnan (y))
498
    return x + y;
499
  if (x == y)
500
    return x;
501
  if (!isfinite (x))
502
    return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
503
 
504
  /* absx = fabsf (x);  */
505
  absx = (x < 0.0) ? -x : x;
506
 
507
  /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals.  */
508
  if (__FLT_DENORM_MIN__ == 0.0f)
509
    denorm_min = __FLT_MIN__;
510
  else
511
    denorm_min = __FLT_DENORM_MIN__;
512
 
513
  if (absx < __FLT_MIN__)
514
    delta = denorm_min;
515
  else
516
    {
517
      float frac;
518
      int exp;
519
 
520
      /* Discard the fraction from x.  */
521
      frac = frexpf (absx, &exp);
522
      delta = scalbnf (0.5f, exp);
523
 
524
      /* Scale x by the epsilon of the representation.  By rights we should
525
         have been able to combine this with scalbnf, but some targets don't
526
         get that correct with denormals.  */
527
      delta *= __FLT_EPSILON__;
528
 
529
      /* If we're going to be reducing the absolute value of X, and doing so
530
         would reduce the exponent of X, then the delta to be applied is
531
         one exponent smaller.  */
532
      if (frac == 0.5f && (y < x) == (x > 0))
533
        delta *= 0.5f;
534
 
535
      /* If that underflows to zero, then we're back to the minimum.  */
536
      if (delta == 0.0f)
537
        delta = denorm_min;
538
    }
539
 
540
  if (y < x)
541
    delta = -delta;
542
 
543
  return x + delta;
544
}
545
#endif
546
 
547
 
548
#if !defined(HAVE_POWF) || defined(HAVE_BROKEN_POWF)
549
#ifndef HAVE_POWF
550
#define HAVE_POWF 1
551
#endif
552
float powf (float x, float y);
553
 
554
float
555
powf (float x, float y)
556
{
557
  return (float) pow (x, y);
558
}
559
#endif
560
 
561
 
562
#ifndef HAVE_ROUND
563
#define HAVE_ROUND 1
564
/* Round to nearest integral value.  If the argument is halfway between two
565
   integral values then round away from zero.  */
566
double round (double x);
567
 
568
double
569
round (double x)
570
{
571
   double t;
572
   if (!isfinite (x))
573
     return (x);
574
 
575
   if (x >= 0.0)
576
    {
577
      t = floor (x);
578
      if (t - x <= -0.5)
579
        t += 1.0;
580
      return (t);
581
    }
582
   else
583
    {
584
      t = floor (-x);
585
      if (t + x <= -0.5)
586
        t += 1.0;
587
      return (-t);
588
    }
589
}
590
#endif
591
 
592
 
593
/* Algorithm by Steven G. Kargl.  */
594
 
595
#if !defined(HAVE_ROUNDL)
596
#define HAVE_ROUNDL 1
597
long double roundl (long double x);
598
 
599
#if defined(HAVE_CEILL)
600
/* Round to nearest integral value.  If the argument is halfway between two
601
   integral values then round away from zero.  */
602
 
603
long double
604
roundl (long double x)
605
{
606
   long double t;
607
   if (!isfinite (x))
608
     return (x);
609
 
610
   if (x >= 0.0)
611
    {
612
      t = ceill (x);
613
      if (t - x > 0.5)
614
        t -= 1.0;
615
      return (t);
616
    }
617
   else
618
    {
619
      t = ceill (-x);
620
      if (t + x > 0.5)
621
        t -= 1.0;
622
      return (-t);
623
    }
624
}
625
#else
626
 
627
/* Poor version of roundl for system that don't have ceill.  */
628
long double
629
roundl (long double x)
630
{
631
  if (x > DBL_MAX || x < -DBL_MAX)
632
    {
633
#ifdef HAVE_NEXTAFTERL
634
      long double prechalf = nextafterl (0.5L, LDBL_MAX);
635
#else
636
      static long double prechalf = 0.5L;
637
#endif
638
      return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf));
639
    }
640
  else
641
    /* Use round().  */
642
    return round ((double) x);
643
}
644
 
645
#endif
646
#endif
647
 
648
#ifndef HAVE_ROUNDF
649
#define HAVE_ROUNDF 1
650
/* Round to nearest integral value.  If the argument is halfway between two
651
   integral values then round away from zero.  */
652
float roundf (float x);
653
 
654
float
655
roundf (float x)
656
{
657
   float t;
658
   if (!isfinite (x))
659
     return (x);
660
 
661
   if (x >= 0.0)
662
    {
663
      t = floorf (x);
664
      if (t - x <= -0.5)
665
        t += 1.0;
666
      return (t);
667
    }
668
   else
669
    {
670
      t = floorf (-x);
671
      if (t + x <= -0.5)
672
        t += 1.0;
673
      return (-t);
674
    }
675
}
676
#endif
677
 
678
 
679
/* lround{f,,l} and llround{f,,l} functions.  */
680
 
681
#if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF)
682
#define HAVE_LROUNDF 1
683
long int lroundf (float x);
684
 
685
long int
686
lroundf (float x)
687
{
688
  return (long int) roundf (x);
689
}
690
#endif
691
 
692
#if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
693
#define HAVE_LROUND 1
694
long int lround (double x);
695
 
696
long int
697
lround (double x)
698
{
699
  return (long int) round (x);
700
}
701
#endif
702
 
703
#if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
704
#define HAVE_LROUNDL 1
705
long int lroundl (long double x);
706
 
707
long int
708
lroundl (long double x)
709
{
710
  return (long long int) roundl (x);
711
}
712
#endif
713
 
714
#if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
715
#define HAVE_LLROUNDF 1
716
long long int llroundf (float x);
717
 
718
long long int
719
llroundf (float x)
720
{
721
  return (long long int) roundf (x);
722
}
723
#endif
724
 
725
#if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
726
#define HAVE_LLROUND 1
727
long long int llround (double x);
728
 
729
long long int
730
llround (double x)
731
{
732
  return (long long int) round (x);
733
}
734
#endif
735
 
736
#if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
737
#define HAVE_LLROUNDL 1
738
long long int llroundl (long double x);
739
 
740
long long int
741
llroundl (long double x)
742
{
743
  return (long long int) roundl (x);
744
}
745
#endif
746
 
747
 
748
#ifndef HAVE_LOG10L
749
#define HAVE_LOG10L 1
750
/* log10 function for long double variables. The version provided here
751
   reduces the argument until it fits into a double, then use log10.  */
752
long double log10l (long double x);
753
 
754
long double
755
log10l (long double x)
756
{
757
#if LDBL_MAX_EXP > DBL_MAX_EXP
758
  if (x > DBL_MAX)
759
    {
760
      double val;
761
      int p2_result = 0;
762
      if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
763
      if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
764
      if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
765
      if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
766
      if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
767
      val = log10 ((double) x);
768
      return (val + p2_result * .30102999566398119521373889472449302L);
769
    }
770
#endif
771
#if LDBL_MIN_EXP < DBL_MIN_EXP
772
  if (x < DBL_MIN)
773
    {
774
      double val;
775
      int p2_result = 0;
776
      if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
777
      if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
778
      if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
779
      if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
780
      if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
781
      val = fabs (log10 ((double) x));
782
      return (- val - p2_result * .30102999566398119521373889472449302L);
783
    }
784
#endif
785
    return log10 (x);
786
}
787
#endif
788
 
789
 
790
#ifndef HAVE_FLOORL
791
#define HAVE_FLOORL 1
792
long double floorl (long double x);
793
 
794
long double
795
floorl (long double x)
796
{
797
  /* Zero, possibly signed.  */
798
  if (x == 0)
799
    return x;
800
 
801
  /* Large magnitude.  */
802
  if (x > DBL_MAX || x < (-DBL_MAX))
803
    return x;
804
 
805
  /* Small positive values.  */
806
  if (x >= 0 && x < DBL_MIN)
807
    return 0;
808
 
809
  /* Small negative values.  */
810
  if (x < 0 && x > (-DBL_MIN))
811
    return -1;
812
 
813
  return floor (x);
814
}
815
#endif
816
 
817
 
818
#ifndef HAVE_FMODL
819
#define HAVE_FMODL 1
820
long double fmodl (long double x, long double y);
821
 
822
long double
823
fmodl (long double x, long double y)
824
{
825
  if (y == 0.0L)
826
    return 0.0L;
827
 
828
  /* Need to check that the result has the same sign as x and magnitude
829
     less than the magnitude of y.  */
830
  return x - floorl (x / y) * y;
831
}
832
#endif
833
 
834
 
835
#if !defined(HAVE_CABSF)
836
#define HAVE_CABSF 1
837
float cabsf (float complex z);
838
 
839
float
840
cabsf (float complex z)
841
{
842
  return hypotf (REALPART (z), IMAGPART (z));
843
}
844
#endif
845
 
846
#if !defined(HAVE_CABS)
847
#define HAVE_CABS 1
848
double cabs (double complex z);
849
 
850
double
851
cabs (double complex z)
852
{
853
  return hypot (REALPART (z), IMAGPART (z));
854
}
855
#endif
856
 
857
#if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
858
#define HAVE_CABSL 1
859
long double cabsl (long double complex z);
860
 
861
long double
862
cabsl (long double complex z)
863
{
864
  return hypotl (REALPART (z), IMAGPART (z));
865
}
866
#endif
867
 
868
 
869
#if !defined(HAVE_CARGF)
870
#define HAVE_CARGF 1
871
float cargf (float complex z);
872
 
873
float
874
cargf (float complex z)
875
{
876
  return atan2f (IMAGPART (z), REALPART (z));
877
}
878
#endif
879
 
880
#if !defined(HAVE_CARG)
881
#define HAVE_CARG 1
882
double carg (double complex z);
883
 
884
double
885
carg (double complex z)
886
{
887
  return atan2 (IMAGPART (z), REALPART (z));
888
}
889
#endif
890
 
891
#if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
892
#define HAVE_CARGL 1
893
long double cargl (long double complex z);
894
 
895
long double
896
cargl (long double complex z)
897
{
898
  return atan2l (IMAGPART (z), REALPART (z));
899
}
900
#endif
901
 
902
 
903
/* exp(z) = exp(a)*(cos(b) + i sin(b))  */
904
#if !defined(HAVE_CEXPF)
905
#define HAVE_CEXPF 1
906
float complex cexpf (float complex z);
907
 
908
float complex
909
cexpf (float complex z)
910
{
911
  float a, b;
912
  float complex v;
913
 
914
  a = REALPART (z);
915
  b = IMAGPART (z);
916
  COMPLEX_ASSIGN (v, cosf (b), sinf (b));
917
  return expf (a) * v;
918
}
919
#endif
920
 
921
#if !defined(HAVE_CEXP)
922
#define HAVE_CEXP 1
923
double complex cexp (double complex z);
924
 
925
double complex
926
cexp (double complex z)
927
{
928
  double a, b;
929
  double complex v;
930
 
931
  a = REALPART (z);
932
  b = IMAGPART (z);
933
  COMPLEX_ASSIGN (v, cos (b), sin (b));
934
  return exp (a) * v;
935
}
936
#endif
937
 
938
#if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
939
#define HAVE_CEXPL 1
940
long double complex cexpl (long double complex z);
941
 
942
long double complex
943
cexpl (long double complex z)
944
{
945
  long double a, b;
946
  long double complex v;
947
 
948
  a = REALPART (z);
949
  b = IMAGPART (z);
950
  COMPLEX_ASSIGN (v, cosl (b), sinl (b));
951
  return expl (a) * v;
952
}
953
#endif
954
 
955
 
956
/* log(z) = log (cabs(z)) + i*carg(z)  */
957
#if !defined(HAVE_CLOGF)
958
#define HAVE_CLOGF 1
959
float complex clogf (float complex z);
960
 
961
float complex
962
clogf (float complex z)
963
{
964
  float complex v;
965
 
966
  COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
967
  return v;
968
}
969
#endif
970
 
971
#if !defined(HAVE_CLOG)
972
#define HAVE_CLOG 1
973
double complex clog (double complex z);
974
 
975
double complex
976
clog (double complex z)
977
{
978
  double complex v;
979
 
980
  COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
981
  return v;
982
}
983
#endif
984
 
985
#if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
986
#define HAVE_CLOGL 1
987
long double complex clogl (long double complex z);
988
 
989
long double complex
990
clogl (long double complex z)
991
{
992
  long double complex v;
993
 
994
  COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
995
  return v;
996
}
997
#endif
998
 
999
 
1000
/* log10(z) = log10 (cabs(z)) + i*carg(z)  */
1001
#if !defined(HAVE_CLOG10F)
1002
#define HAVE_CLOG10F 1
1003
float complex clog10f (float complex z);
1004
 
1005
float complex
1006
clog10f (float complex z)
1007
{
1008
  float complex v;
1009
 
1010
  COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
1011
  return v;
1012
}
1013
#endif
1014
 
1015
#if !defined(HAVE_CLOG10)
1016
#define HAVE_CLOG10 1
1017
double complex clog10 (double complex z);
1018
 
1019
double complex
1020
clog10 (double complex z)
1021
{
1022
  double complex v;
1023
 
1024
  COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
1025
  return v;
1026
}
1027
#endif
1028
 
1029
#if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
1030
#define HAVE_CLOG10L 1
1031
long double complex clog10l (long double complex z);
1032
 
1033
long double complex
1034
clog10l (long double complex z)
1035
{
1036
  long double complex v;
1037
 
1038
  COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
1039
  return v;
1040
}
1041
#endif
1042
 
1043
 
1044
/* pow(base, power) = cexp (power * clog (base))  */
1045
#if !defined(HAVE_CPOWF)
1046
#define HAVE_CPOWF 1
1047
float complex cpowf (float complex base, float complex power);
1048
 
1049
float complex
1050
cpowf (float complex base, float complex power)
1051
{
1052
  return cexpf (power * clogf (base));
1053
}
1054
#endif
1055
 
1056
#if !defined(HAVE_CPOW)
1057
#define HAVE_CPOW 1
1058
double complex cpow (double complex base, double complex power);
1059
 
1060
double complex
1061
cpow (double complex base, double complex power)
1062
{
1063
  return cexp (power * clog (base));
1064
}
1065
#endif
1066
 
1067
#if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
1068
#define HAVE_CPOWL 1
1069
long double complex cpowl (long double complex base, long double complex power);
1070
 
1071
long double complex
1072
cpowl (long double complex base, long double complex power)
1073
{
1074
  return cexpl (power * clogl (base));
1075
}
1076
#endif
1077
 
1078
 
1079
/* sqrt(z).  Algorithm pulled from glibc.  */
1080
#if !defined(HAVE_CSQRTF)
1081
#define HAVE_CSQRTF 1
1082
float complex csqrtf (float complex z);
1083
 
1084
float complex
1085
csqrtf (float complex z)
1086
{
1087
  float re, im;
1088
  float complex v;
1089
 
1090
  re = REALPART (z);
1091
  im = IMAGPART (z);
1092
  if (im == 0)
1093
    {
1094
      if (re < 0)
1095
        {
1096
          COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
1097
        }
1098
      else
1099
        {
1100
          COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
1101
        }
1102
    }
1103
  else if (re == 0)
1104
    {
1105
      float r;
1106
 
1107
      r = sqrtf (0.5 * fabsf (im));
1108
 
1109
      COMPLEX_ASSIGN (v, r, copysignf (r, im));
1110
    }
1111
  else
1112
    {
1113
      float d, r, s;
1114
 
1115
      d = hypotf (re, im);
1116
      /* Use the identity   2  Re res  Im res = Im x
1117
         to avoid cancellation error in  d +/- Re x.  */
1118
      if (re > 0)
1119
        {
1120
          r = sqrtf (0.5 * d + 0.5 * re);
1121
          s = (0.5 * im) / r;
1122
        }
1123
      else
1124
        {
1125
          s = sqrtf (0.5 * d - 0.5 * re);
1126
          r = fabsf ((0.5 * im) / s);
1127
        }
1128
 
1129
      COMPLEX_ASSIGN (v, r, copysignf (s, im));
1130
    }
1131
  return v;
1132
}
1133
#endif
1134
 
1135
#if !defined(HAVE_CSQRT)
1136
#define HAVE_CSQRT 1
1137
double complex csqrt (double complex z);
1138
 
1139
double complex
1140
csqrt (double complex z)
1141
{
1142
  double re, im;
1143
  double complex v;
1144
 
1145
  re = REALPART (z);
1146
  im = IMAGPART (z);
1147
  if (im == 0)
1148
    {
1149
      if (re < 0)
1150
        {
1151
          COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
1152
        }
1153
      else
1154
        {
1155
          COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
1156
        }
1157
    }
1158
  else if (re == 0)
1159
    {
1160
      double r;
1161
 
1162
      r = sqrt (0.5 * fabs (im));
1163
 
1164
      COMPLEX_ASSIGN (v, r, copysign (r, im));
1165
    }
1166
  else
1167
    {
1168
      double d, r, s;
1169
 
1170
      d = hypot (re, im);
1171
      /* Use the identity   2  Re res  Im res = Im x
1172
         to avoid cancellation error in  d +/- Re x.  */
1173
      if (re > 0)
1174
        {
1175
          r = sqrt (0.5 * d + 0.5 * re);
1176
          s = (0.5 * im) / r;
1177
        }
1178
      else
1179
        {
1180
          s = sqrt (0.5 * d - 0.5 * re);
1181
          r = fabs ((0.5 * im) / s);
1182
        }
1183
 
1184
      COMPLEX_ASSIGN (v, r, copysign (s, im));
1185
    }
1186
  return v;
1187
}
1188
#endif
1189
 
1190
#if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
1191
#define HAVE_CSQRTL 1
1192
long double complex csqrtl (long double complex z);
1193
 
1194
long double complex
1195
csqrtl (long double complex z)
1196
{
1197
  long double re, im;
1198
  long double complex v;
1199
 
1200
  re = REALPART (z);
1201
  im = IMAGPART (z);
1202
  if (im == 0)
1203
    {
1204
      if (re < 0)
1205
        {
1206
          COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
1207
        }
1208
      else
1209
        {
1210
          COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
1211
        }
1212
    }
1213
  else if (re == 0)
1214
    {
1215
      long double r;
1216
 
1217
      r = sqrtl (0.5 * fabsl (im));
1218
 
1219
      COMPLEX_ASSIGN (v, copysignl (r, im), r);
1220
    }
1221
  else
1222
    {
1223
      long double d, r, s;
1224
 
1225
      d = hypotl (re, im);
1226
      /* Use the identity   2  Re res  Im res = Im x
1227
         to avoid cancellation error in  d +/- Re x.  */
1228
      if (re > 0)
1229
        {
1230
          r = sqrtl (0.5 * d + 0.5 * re);
1231
          s = (0.5 * im) / r;
1232
        }
1233
      else
1234
        {
1235
          s = sqrtl (0.5 * d - 0.5 * re);
1236
          r = fabsl ((0.5 * im) / s);
1237
        }
1238
 
1239
      COMPLEX_ASSIGN (v, r, copysignl (s, im));
1240
    }
1241
  return v;
1242
}
1243
#endif
1244
 
1245
 
1246
/* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b)  */
1247
#if !defined(HAVE_CSINHF)
1248
#define HAVE_CSINHF 1
1249
float complex csinhf (float complex a);
1250
 
1251
float complex
1252
csinhf (float complex a)
1253
{
1254
  float r, i;
1255
  float complex v;
1256
 
1257
  r = REALPART (a);
1258
  i = IMAGPART (a);
1259
  COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
1260
  return v;
1261
}
1262
#endif
1263
 
1264
#if !defined(HAVE_CSINH)
1265
#define HAVE_CSINH 1
1266
double complex csinh (double complex a);
1267
 
1268
double complex
1269
csinh (double complex a)
1270
{
1271
  double r, i;
1272
  double complex v;
1273
 
1274
  r = REALPART (a);
1275
  i = IMAGPART (a);
1276
  COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
1277
  return v;
1278
}
1279
#endif
1280
 
1281
#if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1282
#define HAVE_CSINHL 1
1283
long double complex csinhl (long double complex a);
1284
 
1285
long double complex
1286
csinhl (long double complex a)
1287
{
1288
  long double r, i;
1289
  long double complex v;
1290
 
1291
  r = REALPART (a);
1292
  i = IMAGPART (a);
1293
  COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1294
  return v;
1295
}
1296
#endif
1297
 
1298
 
1299
/* cosh(a + i b) = cosh(a) cos(b) + i sinh(a) sin(b)  */
1300
#if !defined(HAVE_CCOSHF)
1301
#define HAVE_CCOSHF 1
1302
float complex ccoshf (float complex a);
1303
 
1304
float complex
1305
ccoshf (float complex a)
1306
{
1307
  float r, i;
1308
  float complex v;
1309
 
1310
  r = REALPART (a);
1311
  i = IMAGPART (a);
1312
  COMPLEX_ASSIGN (v, coshf (r) * cosf (i), sinhf (r) * sinf (i));
1313
  return v;
1314
}
1315
#endif
1316
 
1317
#if !defined(HAVE_CCOSH)
1318
#define HAVE_CCOSH 1
1319
double complex ccosh (double complex a);
1320
 
1321
double complex
1322
ccosh (double complex a)
1323
{
1324
  double r, i;
1325
  double complex v;
1326
 
1327
  r = REALPART (a);
1328
  i = IMAGPART (a);
1329
  COMPLEX_ASSIGN (v, cosh (r) * cos (i),  sinh (r) * sin (i));
1330
  return v;
1331
}
1332
#endif
1333
 
1334
#if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1335
#define HAVE_CCOSHL 1
1336
long double complex ccoshl (long double complex a);
1337
 
1338
long double complex
1339
ccoshl (long double complex a)
1340
{
1341
  long double r, i;
1342
  long double complex v;
1343
 
1344
  r = REALPART (a);
1345
  i = IMAGPART (a);
1346
  COMPLEX_ASSIGN (v, coshl (r) * cosl (i), sinhl (r) * sinl (i));
1347
  return v;
1348
}
1349
#endif
1350
 
1351
 
1352
/* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 + i tanh(a) tan(b))  */
1353
#if !defined(HAVE_CTANHF)
1354
#define HAVE_CTANHF 1
1355
float complex ctanhf (float complex a);
1356
 
1357
float complex
1358
ctanhf (float complex a)
1359
{
1360
  float rt, it;
1361
  float complex n, d;
1362
 
1363
  rt = tanhf (REALPART (a));
1364
  it = tanf (IMAGPART (a));
1365
  COMPLEX_ASSIGN (n, rt, it);
1366
  COMPLEX_ASSIGN (d, 1, rt * it);
1367
 
1368
  return n / d;
1369
}
1370
#endif
1371
 
1372
#if !defined(HAVE_CTANH)
1373
#define HAVE_CTANH 1
1374
double complex ctanh (double complex a);
1375
double complex
1376
ctanh (double complex a)
1377
{
1378
  double rt, it;
1379
  double complex n, d;
1380
 
1381
  rt = tanh (REALPART (a));
1382
  it = tan (IMAGPART (a));
1383
  COMPLEX_ASSIGN (n, rt, it);
1384
  COMPLEX_ASSIGN (d, 1, rt * it);
1385
 
1386
  return n / d;
1387
}
1388
#endif
1389
 
1390
#if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1391
#define HAVE_CTANHL 1
1392
long double complex ctanhl (long double complex a);
1393
 
1394
long double complex
1395
ctanhl (long double complex a)
1396
{
1397
  long double rt, it;
1398
  long double complex n, d;
1399
 
1400
  rt = tanhl (REALPART (a));
1401
  it = tanl (IMAGPART (a));
1402
  COMPLEX_ASSIGN (n, rt, it);
1403
  COMPLEX_ASSIGN (d, 1, rt * it);
1404
 
1405
  return n / d;
1406
}
1407
#endif
1408
 
1409
 
1410
/* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b)  */
1411
#if !defined(HAVE_CSINF)
1412
#define HAVE_CSINF 1
1413
float complex csinf (float complex a);
1414
 
1415
float complex
1416
csinf (float complex a)
1417
{
1418
  float r, i;
1419
  float complex v;
1420
 
1421
  r = REALPART (a);
1422
  i = IMAGPART (a);
1423
  COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1424
  return v;
1425
}
1426
#endif
1427
 
1428
#if !defined(HAVE_CSIN)
1429
#define HAVE_CSIN 1
1430
double complex csin (double complex a);
1431
 
1432
double complex
1433
csin (double complex a)
1434
{
1435
  double r, i;
1436
  double complex v;
1437
 
1438
  r = REALPART (a);
1439
  i = IMAGPART (a);
1440
  COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1441
  return v;
1442
}
1443
#endif
1444
 
1445
#if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1446
#define HAVE_CSINL 1
1447
long double complex csinl (long double complex a);
1448
 
1449
long double complex
1450
csinl (long double complex a)
1451
{
1452
  long double r, i;
1453
  long double complex v;
1454
 
1455
  r = REALPART (a);
1456
  i = IMAGPART (a);
1457
  COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1458
  return v;
1459
}
1460
#endif
1461
 
1462
 
1463
/* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b)  */
1464
#if !defined(HAVE_CCOSF)
1465
#define HAVE_CCOSF 1
1466
float complex ccosf (float complex a);
1467
 
1468
float complex
1469
ccosf (float complex a)
1470
{
1471
  float r, i;
1472
  float complex v;
1473
 
1474
  r = REALPART (a);
1475
  i = IMAGPART (a);
1476
  COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1477
  return v;
1478
}
1479
#endif
1480
 
1481
#if !defined(HAVE_CCOS)
1482
#define HAVE_CCOS 1
1483
double complex ccos (double complex a);
1484
 
1485
double complex
1486
ccos (double complex a)
1487
{
1488
  double r, i;
1489
  double complex v;
1490
 
1491
  r = REALPART (a);
1492
  i = IMAGPART (a);
1493
  COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1494
  return v;
1495
}
1496
#endif
1497
 
1498
#if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1499
#define HAVE_CCOSL 1
1500
long double complex ccosl (long double complex a);
1501
 
1502
long double complex
1503
ccosl (long double complex a)
1504
{
1505
  long double r, i;
1506
  long double complex v;
1507
 
1508
  r = REALPART (a);
1509
  i = IMAGPART (a);
1510
  COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1511
  return v;
1512
}
1513
#endif
1514
 
1515
 
1516
/* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b))  */
1517
#if !defined(HAVE_CTANF)
1518
#define HAVE_CTANF 1
1519
float complex ctanf (float complex a);
1520
 
1521
float complex
1522
ctanf (float complex a)
1523
{
1524
  float rt, it;
1525
  float complex n, d;
1526
 
1527
  rt = tanf (REALPART (a));
1528
  it = tanhf (IMAGPART (a));
1529
  COMPLEX_ASSIGN (n, rt, it);
1530
  COMPLEX_ASSIGN (d, 1, - (rt * it));
1531
 
1532
  return n / d;
1533
}
1534
#endif
1535
 
1536
#if !defined(HAVE_CTAN)
1537
#define HAVE_CTAN 1
1538
double complex ctan (double complex a);
1539
 
1540
double complex
1541
ctan (double complex a)
1542
{
1543
  double rt, it;
1544
  double complex n, d;
1545
 
1546
  rt = tan (REALPART (a));
1547
  it = tanh (IMAGPART (a));
1548
  COMPLEX_ASSIGN (n, rt, it);
1549
  COMPLEX_ASSIGN (d, 1, - (rt * it));
1550
 
1551
  return n / d;
1552
}
1553
#endif
1554
 
1555
#if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1556
#define HAVE_CTANL 1
1557
long double complex ctanl (long double complex a);
1558
 
1559
long double complex
1560
ctanl (long double complex a)
1561
{
1562
  long double rt, it;
1563
  long double complex n, d;
1564
 
1565
  rt = tanl (REALPART (a));
1566
  it = tanhl (IMAGPART (a));
1567
  COMPLEX_ASSIGN (n, rt, it);
1568
  COMPLEX_ASSIGN (d, 1, - (rt * it));
1569
 
1570
  return n / d;
1571
}
1572
#endif
1573
 
1574
 
1575
/* Complex ASIN.  Returns wrongly NaN for infinite arguments.
1576
   Algorithm taken from Abramowitz & Stegun.  */
1577
 
1578
#if !defined(HAVE_CASINF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1579
#define HAVE_CASINF 1
1580
complex float casinf (complex float z);
1581
 
1582
complex float
1583
casinf (complex float z)
1584
{
1585
  return -I*clogf (I*z + csqrtf (1.0f-z*z));
1586
}
1587
#endif
1588
 
1589
 
1590
#if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1591
#define HAVE_CASIN 1
1592
complex double casin (complex double z);
1593
 
1594
complex double
1595
casin (complex double z)
1596
{
1597
  return -I*clog (I*z + csqrt (1.0-z*z));
1598
}
1599
#endif
1600
 
1601
 
1602
#if !defined(HAVE_CASINL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1603
#define HAVE_CASINL 1
1604
complex long double casinl (complex long double z);
1605
 
1606
complex long double
1607
casinl (complex long double z)
1608
{
1609
  return -I*clogl (I*z + csqrtl (1.0L-z*z));
1610
}
1611
#endif
1612
 
1613
 
1614
/* Complex ACOS.  Returns wrongly NaN for infinite arguments.
1615
   Algorithm taken from Abramowitz & Stegun.  */
1616
 
1617
#if !defined(HAVE_CACOSF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1618
#define HAVE_CACOSF 1
1619
complex float cacosf (complex float z);
1620
 
1621
complex float
1622
cacosf (complex float z)
1623
{
1624
  return -I*clogf (z + I*csqrtf (1.0f-z*z));
1625
}
1626
#endif
1627
 
1628
 
1629
#if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1630
#define HAVE_CACOS 1
1631
complex double cacos (complex double z);
1632
 
1633
complex double
1634
cacos (complex double z)
1635
{
1636
  return -I*clog (z + I*csqrt (1.0-z*z));
1637
}
1638
#endif
1639
 
1640
 
1641
#if !defined(HAVE_CACOSL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1642
#define HAVE_CACOSL 1
1643
complex long double cacosl (complex long double z);
1644
 
1645
complex long double
1646
cacosl (complex long double z)
1647
{
1648
  return -I*clogl (z + I*csqrtl (1.0L-z*z));
1649
}
1650
#endif
1651
 
1652
 
1653
/* Complex ATAN.  Returns wrongly NaN for infinite arguments.
1654
   Algorithm taken from Abramowitz & Stegun.  */
1655
 
1656
#if !defined(HAVE_CATANF) && defined(HAVE_CLOGF)
1657
#define HAVE_CACOSF 1
1658
complex float catanf (complex float z);
1659
 
1660
complex float
1661
catanf (complex float z)
1662
{
1663
  return I*clogf ((I+z)/(I-z))/2.0f;
1664
}
1665
#endif
1666
 
1667
 
1668
#if !defined(HAVE_CATAN) && defined(HAVE_CLOG)
1669
#define HAVE_CACOS 1
1670
complex double catan (complex double z);
1671
 
1672
complex double
1673
catan (complex double z)
1674
{
1675
  return I*clog ((I+z)/(I-z))/2.0;
1676
}
1677
#endif
1678
 
1679
 
1680
#if !defined(HAVE_CATANL) && defined(HAVE_CLOGL)
1681
#define HAVE_CACOSL 1
1682
complex long double catanl (complex long double z);
1683
 
1684
complex long double
1685
catanl (complex long double z)
1686
{
1687
  return I*clogl ((I+z)/(I-z))/2.0L;
1688
}
1689
#endif
1690
 
1691
 
1692
/* Complex ASINH.  Returns wrongly NaN for infinite arguments.
1693
   Algorithm taken from Abramowitz & Stegun.  */
1694
 
1695
#if !defined(HAVE_CASINHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1696
#define HAVE_CASINHF 1
1697
complex float casinhf (complex float z);
1698
 
1699
complex float
1700
casinhf (complex float z)
1701
{
1702
  return clogf (z + csqrtf (z*z+1.0f));
1703
}
1704
#endif
1705
 
1706
 
1707
#if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1708
#define HAVE_CASINH 1
1709
complex double casinh (complex double z);
1710
 
1711
complex double
1712
casinh (complex double z)
1713
{
1714
  return clog (z + csqrt (z*z+1.0));
1715
}
1716
#endif
1717
 
1718
 
1719
#if !defined(HAVE_CASINHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1720
#define HAVE_CASINHL 1
1721
complex long double casinhl (complex long double z);
1722
 
1723
complex long double
1724
casinhl (complex long double z)
1725
{
1726
  return clogl (z + csqrtl (z*z+1.0L));
1727
}
1728
#endif
1729
 
1730
 
1731
/* Complex ACOSH.  Returns wrongly NaN for infinite arguments.
1732
   Algorithm taken from Abramowitz & Stegun.  */
1733
 
1734
#if !defined(HAVE_CACOSHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1735
#define HAVE_CACOSHF 1
1736
complex float cacoshf (complex float z);
1737
 
1738
complex float
1739
cacoshf (complex float z)
1740
{
1741
  return clogf (z + csqrtf (z-1.0f) * csqrtf (z+1.0f));
1742
}
1743
#endif
1744
 
1745
 
1746
#if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1747
#define HAVE_CACOSH 1
1748
complex double cacosh (complex double z);
1749
 
1750
complex double
1751
cacosh (complex double z)
1752
{
1753
  return clog (z + csqrt (z-1.0) * csqrt (z+1.0));
1754
}
1755
#endif
1756
 
1757
 
1758
#if !defined(HAVE_CACOSHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1759
#define HAVE_CACOSHL 1
1760
complex long double cacoshl (complex long double z);
1761
 
1762
complex long double
1763
cacoshl (complex long double z)
1764
{
1765
  return clogl (z + csqrtl (z-1.0L) * csqrtl (z+1.0L));
1766
}
1767
#endif
1768
 
1769
 
1770
/* Complex ATANH.  Returns wrongly NaN for infinite arguments.
1771
   Algorithm taken from Abramowitz & Stegun.  */
1772
 
1773
#if !defined(HAVE_CATANHF) && defined(HAVE_CLOGF)
1774
#define HAVE_CATANHF 1
1775
complex float catanhf (complex float z);
1776
 
1777
complex float
1778
catanhf (complex float z)
1779
{
1780
  return clogf ((1.0f+z)/(1.0f-z))/2.0f;
1781
}
1782
#endif
1783
 
1784
 
1785
#if !defined(HAVE_CATANH) && defined(HAVE_CLOG)
1786
#define HAVE_CATANH 1
1787
complex double catanh (complex double z);
1788
 
1789
complex double
1790
catanh (complex double z)
1791
{
1792
  return clog ((1.0+z)/(1.0-z))/2.0;
1793
}
1794
#endif
1795
 
1796
#if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL)
1797
#define HAVE_CATANHL 1
1798
complex long double catanhl (complex long double z);
1799
 
1800
complex long double
1801
catanhl (complex long double z)
1802
{
1803
  return clogl ((1.0L+z)/(1.0L-z))/2.0L;
1804
}
1805
#endif
1806
 
1807
 
1808
#if !defined(HAVE_TGAMMA)
1809
#define HAVE_TGAMMA 1
1810
double tgamma (double);
1811
 
1812
/* Fallback tgamma() function. Uses the algorithm from
1813
   http://www.netlib.org/specfun/gamma and references therein.  */
1814
 
1815
#undef SQRTPI
1816
#define SQRTPI 0.9189385332046727417803297
1817
 
1818
#undef PI
1819
#define PI 3.1415926535897932384626434
1820
 
1821
double
1822
tgamma (double x)
1823
{
1824
  int i, n, parity;
1825
  double fact, res, sum, xden, xnum, y, y1, ysq, z;
1826
 
1827
  static double p[8] = {
1828
    -1.71618513886549492533811e0,  2.47656508055759199108314e1,
1829
    -3.79804256470945635097577e2,  6.29331155312818442661052e2,
1830
     8.66966202790413211295064e2, -3.14512729688483675254357e4,
1831
    -3.61444134186911729807069e4,  6.64561438202405440627855e4 };
1832
 
1833
  static double q[8] = {
1834
    -3.08402300119738975254353e1,  3.15350626979604161529144e2,
1835
    -1.01515636749021914166146e3, -3.10777167157231109440444e3,
1836
     2.25381184209801510330112e4,  4.75584627752788110767815e3,
1837
    -1.34659959864969306392456e5, -1.15132259675553483497211e5 };
1838
 
1839
  static double c[7] = {             -1.910444077728e-03,
1840
     8.4171387781295e-04,            -5.952379913043012e-04,
1841
     7.93650793500350248e-04,        -2.777777777777681622553e-03,
1842
     8.333333333333333331554247e-02,  5.7083835261e-03 };
1843
 
1844
  static const double xminin = 2.23e-308;
1845
  static const double xbig = 171.624;
1846
  static const double xnan = __builtin_nan ("0x0"), xinf = __builtin_inf ();
1847
  static double eps = 0;
1848
 
1849
  if (eps == 0)
1850
    eps = nextafter (1., 2.) - 1.;
1851
 
1852
  parity = 0;
1853
  fact = 1;
1854
  n = 0;
1855
  y = x;
1856
 
1857
  if (isnan (x))
1858
    return x;
1859
 
1860
  if (y <= 0)
1861
    {
1862
      y = -x;
1863
      y1 = trunc (y);
1864
      res = y - y1;
1865
 
1866
      if (res != 0)
1867
        {
1868
          if (y1 != trunc (y1*0.5l)*2)
1869
            parity = 1;
1870
          fact = -PI / sin (PI*res);
1871
          y = y + 1;
1872
        }
1873
      else
1874
        return x == 0 ? copysign (xinf, x) : xnan;
1875
    }
1876
 
1877
  if (y < eps)
1878
    {
1879
      if (y >= xminin)
1880
        res = 1 / y;
1881
      else
1882
        return xinf;
1883
    }
1884
  else if (y < 13)
1885
    {
1886
      y1 = y;
1887
      if (y < 1)
1888
        {
1889
          z = y;
1890
          y = y + 1;
1891
        }
1892
      else
1893
        {
1894
          n = (int)y - 1;
1895
          y = y - n;
1896
          z = y - 1;
1897
        }
1898
 
1899
      xnum = 0;
1900
      xden = 1;
1901
      for (i = 0; i < 8; i++)
1902
        {
1903
          xnum = (xnum + p[i]) * z;
1904
          xden = xden * z + q[i];
1905
        }
1906
 
1907
      res = xnum / xden + 1;
1908
 
1909
      if (y1 < y)
1910
        res = res / y1;
1911
      else if (y1 > y)
1912
        for (i = 1; i <= n; i++)
1913
          {
1914
            res = res * y;
1915
            y = y + 1;
1916
          }
1917
    }
1918
  else
1919
    {
1920
      if (y < xbig)
1921
        {
1922
          ysq = y * y;
1923
          sum = c[6];
1924
          for (i = 0; i < 6; i++)
1925
            sum = sum / ysq + c[i];
1926
 
1927
          sum = sum/y - y + SQRTPI;
1928
          sum = sum + (y - 0.5) * log (y);
1929
          res = exp (sum);
1930
        }
1931
      else
1932
        return x < 0 ? xnan : xinf;
1933
    }
1934
 
1935
  if (parity)
1936
    res = -res;
1937
  if (fact != 1)
1938
    res = fact / res;
1939
 
1940
  return res;
1941
}
1942
#endif
1943
 
1944
 
1945
 
1946
#if !defined(HAVE_LGAMMA)
1947
#define HAVE_LGAMMA 1
1948
double lgamma (double);
1949
 
1950
/* Fallback lgamma() function. Uses the algorithm from
1951
   http://www.netlib.org/specfun/algama and references therein,
1952
   except for negative arguments (where netlib would return +Inf)
1953
   where we use the following identity:
1954
       lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
1955
 */
1956
 
1957
double
1958
lgamma (double y)
1959
{
1960
 
1961
#undef SQRTPI
1962
#define SQRTPI 0.9189385332046727417803297
1963
 
1964
#undef PI
1965
#define PI 3.1415926535897932384626434
1966
 
1967
#define PNT68  0.6796875
1968
#define D1    -0.5772156649015328605195174
1969
#define D2     0.4227843350984671393993777
1970
#define D4     1.791759469228055000094023
1971
 
1972
  static double p1[8] = {
1973
              4.945235359296727046734888e0, 2.018112620856775083915565e2,
1974
              2.290838373831346393026739e3, 1.131967205903380828685045e4,
1975
              2.855724635671635335736389e4, 3.848496228443793359990269e4,
1976
              2.637748787624195437963534e4, 7.225813979700288197698961e3 };
1977
  static double q1[8] = {
1978
              6.748212550303777196073036e1,  1.113332393857199323513008e3,
1979
              7.738757056935398733233834e3,  2.763987074403340708898585e4,
1980
              5.499310206226157329794414e4,  6.161122180066002127833352e4,
1981
              3.635127591501940507276287e4,  8.785536302431013170870835e3 };
1982
  static double p2[8] = {
1983
              4.974607845568932035012064e0,  5.424138599891070494101986e2,
1984
              1.550693864978364947665077e4,  1.847932904445632425417223e5,
1985
              1.088204769468828767498470e6,  3.338152967987029735917223e6,
1986
              5.106661678927352456275255e6,  3.074109054850539556250927e6 };
1987
  static double q2[8] = {
1988
              1.830328399370592604055942e2,  7.765049321445005871323047e3,
1989
              1.331903827966074194402448e5,  1.136705821321969608938755e6,
1990
              5.267964117437946917577538e6,  1.346701454311101692290052e7,
1991
              1.782736530353274213975932e7,  9.533095591844353613395747e6 };
1992
  static double p4[8] = {
1993
              1.474502166059939948905062e4,  2.426813369486704502836312e6,
1994
              1.214755574045093227939592e8,  2.663432449630976949898078e9,
1995
              2.940378956634553899906876e10, 1.702665737765398868392998e11,
1996
              4.926125793377430887588120e11, 5.606251856223951465078242e11 };
1997
  static double q4[8] = {
1998
              2.690530175870899333379843e3,  6.393885654300092398984238e5,
1999
              4.135599930241388052042842e7,  1.120872109616147941376570e9,
2000
              1.488613728678813811542398e10, 1.016803586272438228077304e11,
2001
              3.417476345507377132798597e11, 4.463158187419713286462081e11 };
2002
  static double  c[7] = {
2003
             -1.910444077728e-03,            8.4171387781295e-04,
2004
             -5.952379913043012e-04,         7.93650793500350248e-04,
2005
             -2.777777777777681622553e-03,   8.333333333333333331554247e-02,
2006
              5.7083835261e-03 };
2007
 
2008
  static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0,
2009
                frtbig = 2.25e76;
2010
 
2011
  int i;
2012
  double corr, res, xden, xm1, xm2, xm4, xnum, ysq;
2013
 
2014
  if (eps == 0)
2015
    eps = __builtin_nextafter (1., 2.) - 1.;
2016
 
2017
  if ((y > 0) && (y <= xbig))
2018
    {
2019
      if (y <= eps)
2020
        res = -log (y);
2021
      else if (y <= 1.5)
2022
        {
2023
          if (y < PNT68)
2024
            {
2025
              corr = -log (y);
2026
              xm1 = y;
2027
            }
2028
          else
2029
            {
2030
              corr = 0;
2031
              xm1 = (y - 0.5) - 0.5;
2032
            }
2033
 
2034
          if ((y <= 0.5) || (y >= PNT68))
2035
            {
2036
              xden = 1;
2037
              xnum = 0;
2038
              for (i = 0; i < 8; i++)
2039
                {
2040
                  xnum = xnum*xm1 + p1[i];
2041
                  xden = xden*xm1 + q1[i];
2042
                }
2043
              res = corr + (xm1 * (D1 + xm1*(xnum/xden)));
2044
            }
2045
          else
2046
            {
2047
              xm2 = (y - 0.5) - 0.5;
2048
              xden = 1;
2049
              xnum = 0;
2050
              for (i = 0; i < 8; i++)
2051
                {
2052
                  xnum = xnum*xm2 + p2[i];
2053
                  xden = xden*xm2 + q2[i];
2054
                }
2055
              res = corr + xm2 * (D2 + xm2*(xnum/xden));
2056
            }
2057
        }
2058
      else if (y <= 4)
2059
        {
2060
          xm2 = y - 2;
2061
          xden = 1;
2062
          xnum = 0;
2063
          for (i = 0; i < 8; i++)
2064
            {
2065
              xnum = xnum*xm2 + p2[i];
2066
              xden = xden*xm2 + q2[i];
2067
            }
2068
          res = xm2 * (D2 + xm2*(xnum/xden));
2069
        }
2070
      else if (y <= 12)
2071
        {
2072
          xm4 = y - 4;
2073
          xden = -1;
2074
          xnum = 0;
2075
          for (i = 0; i < 8; i++)
2076
            {
2077
              xnum = xnum*xm4 + p4[i];
2078
              xden = xden*xm4 + q4[i];
2079
            }
2080
          res = D4 + xm4*(xnum/xden);
2081
        }
2082
      else
2083
        {
2084
          res = 0;
2085
          if (y <= frtbig)
2086
            {
2087
              res = c[6];
2088
              ysq = y * y;
2089
              for (i = 0; i < 6; i++)
2090
                res = res / ysq + c[i];
2091
            }
2092
          res = res/y;
2093
          corr = log (y);
2094
          res = res + SQRTPI - 0.5*corr;
2095
          res = res + y*(corr-1);
2096
        }
2097
    }
2098
  else if (y < 0 && __builtin_floor (y) != y)
2099
    {
2100
      /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
2101
         For abs(y) very close to zero, we use a series expansion to
2102
         the first order in y to avoid overflow.  */
2103
      if (y > -1.e-100)
2104
        res = -2 * log (fabs (y)) - lgamma (-y);
2105
      else
2106
        res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y);
2107
    }
2108
  else
2109
    res = xinf;
2110
 
2111
  return res;
2112
}
2113
#endif
2114
 
2115
 
2116
#if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF)
2117
#define HAVE_TGAMMAF 1
2118
float tgammaf (float);
2119
 
2120
float
2121
tgammaf (float x)
2122
{
2123
  return (float) tgamma ((double) x);
2124
}
2125
#endif
2126
 
2127
#if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF)
2128
#define HAVE_LGAMMAF 1
2129
float lgammaf (float);
2130
 
2131
float
2132
lgammaf (float x)
2133
{
2134
  return (float) lgamma ((double) x);
2135
}
2136
#endif

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.