overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22//
23// AuxiliaryDD.c
24//
25// Double-double Function Library
26// Copyright: � 1995-96 by Apple Computer, Inc., all rights reserved
27//
28
29#include "math.h"
30#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
31#include "DD.h"
32#include "fp_private.h"
33#include "limits.h"
34
35extern long double sqrtl(long double);
36
37#define REM_NAN "9"
38#define kSignMask 0x80000000u
39#define kExpMask 0x7ff00000u
40
41static const Double kTZ = {{0x0,0x1}};
42static const Double kUP = {{0x0,0x2}};
43static const Double kDN = {{0x0,0x3}};
44static const float kTwoTo52 = 4503599627370496.0f; // 0x1.0p+52
45
46long double logbl( long double x )
47{
48 DblDbl xDD;
49 double fpenv, sum;
50
51 xDD.f = x;
52
53 FEGETENVD(fpenv);
54 FESETENVD(kTZ.f); // round toward zero
55 sum = xDD.d[0] + xDD.d[1]; // sum head and tail
56 FESETENVD(fpenv);
57
58 return logb(sum);
59}
60
61int ilogbl( long double x )
62{
63 DblDbl xDD;
64 double fpenv, sum;
65
66 xDD.f = x;
67
68 FEGETENVD(fpenv);
69 FESETENVD(kTZ.f); // round toward zero
70 sum = xDD.d[0] + xDD.d[1]; // sum head and tail
71 FESETENVD(fpenv);
72
73 return ilogb(sum);
74}
75
76long double scalbnl( long double x, int iscale )
77{
78 DblDbl xDD;
79 double fpenv;
80 double Z,z;
81 uint32_t hibits;
82
83 FEGETENVD(fpenv);
84 FESETENVD(0.0);
85 xDD.f = x;
86 xDD.d[0] = scalbn(xDD.d[0],iscale); // scale head
87 hibits = (xDD.i[0] >> 20) & 0x7ffu; // biased exp of scaled head
88 if ((hibits != 0u) && (hibits != 0x7ffu)) { // normal case
89 z = scalbn(xDD.d[1],iscale); // scale tail
90 FESETENVD(kTZ.f); // round toward zero
91 Z = xDD.d[0] + z; // get in canonical form
92 xDD.d[1] = (xDD.d[0] - Z) + z;
93 xDD.d[0] = Z;
94 FESETENVD(fpenv);
95 return (xDD.f); // return normal result
96 }
97
98 // head of result is infinite, zero, subnormal, or NaN
99 Z = xDD.d[0];
100 if ((Z != Z) || (Z == 0.0)) // NaN or zero result
101 xDD.d[1] = Z;
102 else if (hibits == 0x7ffu) // infinite result
103 xDD.d[1] = 0.0;
104 else // subnormal result
105 xDD.d[1] = scalbn(xDD.d[1],iscale);
106
107 FESETENVD(fpenv);
108 return (xDD.f); // return result
109}
110
111long double scalblnl( long double x, long int n )
112{
113 int m;
114
115 // Clip n
116 if (unlikely(n > 2097))
117 m = 2098;
118 else if (unlikely(n < -2098))
119 m = -2099;
120 else
121 m = (int) n;
122
123 return scalbnl(x, m);
124}
125
126// These work just as well for the LP64 ABI
127long int lrintl ( long double x )
128{
129 long double t;
130 long int result;
131 double fenv;
132
133 if (unlikely(x != x))
134 {
135 feraiseexcept(FE_INVALID);
136 return LONG_MIN;
137 }
138
139 FEGETENVD(fenv);
140 t = rintl ( x );
141 FESETENVD(fenv);
142
143 if ( t < (long double)LONG_MIN )
144 {
145 feraiseexcept(FE_INVALID);
146 result = LONG_MIN;
147 }
148 else if ( t > (long double)LONG_MAX )
149 {
150 feraiseexcept(FE_INVALID);
151 result = LONG_MAX;
152 }
153 else if (t != x)
154 {
155 feraiseexcept(FE_INEXACT);
156 result = (long int) t;
157 }
158 else
159 {
160 result = (long int) t;
161 }
162
163 return result;
164}
165
166long double hypotl( long double x, long double y )
167{
168 DblDbl xDD, yDD;
169 double fpenv;
170 long double ldtemp;
171 uint32_t expx, expy;
172 const Float kINF = {0x7f800000};
173
174 FEGETENVD(fpenv);
175 FESETENVD(0.0);
176
177 xDD.f = x;
178 yDD.f = y;
179
180 // calculate absolute values of x and y
181 if (xDD.i[0] >= kSignMask) {
182 xDD.i[0] ^= kSignMask;
183 xDD.i[2] ^= kSignMask;
184 }
185 if (yDD.i[0] >= kSignMask) {
186 yDD.i[0] ^= kSignMask;
187 yDD.i[2] ^= kSignMask;
188 }
189
190 expx = (xDD.i[0] >> 20) & 0x7ffu;
191 expy = (yDD.i[0] >> 20) & 0x7ffu;
192
193 // NaN or INF arg(s)
194 if ((expx == 0x7ffu) || (expy == 0x7ffu)) {
195 if (xDD.d[0] == kINF.f) { // propagate INF
196 FESETENVD(fpenv);
197 return (xDD.f);
198 }
199 else if (yDD.d[0] == kINF.f) {
200 FESETENVD(fpenv);
201 return (yDD.f);
202 }
203 else {
204 ldtemp = x + y;
205 FESETENVD(fpenv);
206 return ldtemp; // propagate NaN
207 }
208 }
209
210 // finite arguments at this point
211 if (yDD.f > xDD.f) { // order arguments
212 ldtemp = yDD.f;
213 yDD.f = xDD.f;
214 xDD.f = ldtemp;
215 }
216 if (yDD.d[0] == 0.0) { // zero argument(s)
217 FESETENVD(fpenv);
218 return (xDD.f);
219 }
220
221 // usual case
222 ldtemp = yDD.f / xDD.f; // (avoids UNDERFLOW)
223 ldtemp *= ldtemp;
224 ldtemp = xDD.f*sqrtl(1.0L + ldtemp);
225
226 FESETENVD(fpenv);
227 return (ldtemp);
228}
229
230
231long double truncl( long double x )
232{
233 DblDbl ldu;
234 double xd, result, fpenv;
235
236 FEGETENVD(fpenv);
237 FESETENVD(0.0);
238
239 ldu.f = x;
240
241 if ((ldu.i[0] & kExpMask) == kExpMask) { // NaN or INF
242 FESETENVD(fpenv);
243 return (x);
244 }
245
246 if (fabs(ldu.d[1]) >= kTwoTo52) { // large integral x
247 FESETENVD(fpenv);
248 return (x);
249 }
250
251 // binary point is within or to the left of x's bits
252 FESETENVD(kTZ.f); // round toward zero
253 xd = ldu.d[0] + ldu.d[1]; // sum head and tail
254
255 if (fabs(xd) < kTwoTo52) { // binary point in head
256 ldu.d[1] = 0.0; // result tail must be zero
257 if (xd >= 0.0) // round head to integral
258 result = (xd + kTwoTo52) - kTwoTo52;
259 else
260 result = (xd - kTwoTo52) + kTwoTo52;
261 if (result == 0.0) { // preserve sign of zero
262 if (ldu.i[0] < kSignMask)
263 result = +0.0;
264 else {
265 result = -0.0;
266 ldu.d[1] = -0.0;
267 }
268 }
269 ldu.d[0] = result;
270 FESETENVD(fpenv); // restore caller's rounding
271 return (ldu.f); // deliver result
272 }
273
274 // binary point in tail or between head and tail
275 if (xd > 0.0) // round tail to integral
276 {
277 FESETENVD(kDN.f);
278 }
279 else
280 {
281 FESETENVD(kUP.f);
282 }
283 if (ldu.d[1] >= 0.0)
284 result = (ldu.d[1] + kTwoTo52) - kTwoTo52;
285 else
286 result = (ldu.d[1] - kTwoTo52) + kTwoTo52;
287
288 // make result "canonical" in to-nearest rounding (preserves value)
289 FESETENVD(0.0);
290
291 xd = ldu.d[0] + result;
292 ldu.d[1] = (ldu.d[0] - xd) + result;
293 ldu.d[0] = xd;
294
295 FESETENVD(fpenv);
296 return (ldu.f); // deliver result
297}
298
299
300long double floorl( long double x )
301{
302 DblDbl ldu;
303 double xd, result, fpenv;
304
305 FEGETENVD(fpenv);
306 FESETENVD(0.0);
307
308 ldu.f = x;
309
310 if ((ldu.i[0] & kExpMask) == kExpMask) { // NaN or INF
311 FESETENVD(fpenv);
312 return (x);
313 }
314
315 if (fabs(ldu.d[1]) >= kTwoTo52) { // large integral x
316 FESETENVD(fpenv);
317 return (x);
318 }
319
320 // binary point is within or to the left of x's bits
321 FESETENVD(kDN.f); // round downward
322
323 if (fabs(ldu.d[0]) < kTwoTo52) { // binary point in head
324 xd = ldu.d[0] + ldu.d[1]; // sum head and tail of x
325 ldu.d[1] = 0.0; // result tail is zero
326 if (ldu.d[0] >= 0.0) // round to integral
327 result = (xd + kTwoTo52) - kTwoTo52;
328 else
329 result = (xd - kTwoTo52) + kTwoTo52;
330 if (result == 0.0) { // preserve sign of zero
331 if (ldu.i[0] < kSignMask)
332 result = 0.0;
333 else {
334 result = -0.0;
335 ldu.d[1] = result;
336 }
337 }
338 ldu.d[0] = result;
339 FESETENVD(fpenv); // restore caller's rounding
340 return (ldu.f); // deliver result
341 }
342
343 // binary point in tail or between head and tail;
344 if (ldu.d[1] >= 0.0) // round to integral
345 result = (ldu.d[1] + kTwoTo52) - kTwoTo52;
346 else
347 result = (ldu.d[1] - kTwoTo52) + kTwoTo52;
348
349 // make result "canonical" in to-nearest rounding (preserves value)
350 FESETENVD(0.0);
351
352 xd = ldu.d[0] + result;
353 ldu.d[1] = (ldu.d[0] - xd) + result;
354 ldu.d[0] = xd;
355
356 FESETENVD(fpenv);
357 return (ldu.f); // deliver result
358}
359
360
361long double ceill( long double x )
362{
363 DblDbl ldu;
364 double xd, result, fpenv;
365
366 FEGETENVD(fpenv);
367 FESETENVD(0.0);
368
369 ldu.f = x;
370
371 if ((ldu.i[0] & kExpMask) == kExpMask) { // NaN or INF
372 FESETENVD(fpenv);
373 return (x);
374 }
375
376 if (fabs(ldu.d[1]) >= kTwoTo52) { // large integral x
377 FESETENVD(fpenv);
378 return (x);
379 }
380
381 // binary point is within or to the left of x's bits
382 FESETENVD(kUP.f); // round upward
383
384 if (fabs(ldu.d[0]) < kTwoTo52) { // binary point in head
385 xd = ldu.d[0] + ldu.d[1]; // sum head and tail of x
386 ldu.d[1] = 0.0; // result tail is zero
387 if (ldu.d[0] >= 0.0) // round to integral
388 result = (xd + kTwoTo52) - kTwoTo52;
389 else
390 result = (xd - kTwoTo52) + kTwoTo52;
391 if (result == 0.0) { // preserve sign of zero
392 if (ldu.i[0] < kSignMask)
393 result = 0.0;
394 else {
395 result = -0.0;
396 ldu.d[1] = result;
397 }
398 }
399 ldu.d[0] = result;
400 FESETENVD(fpenv);
401 // restore caller's rounding
402 return (ldu.f); // deliver result
403 }
404
405 // binary point in tail or between head and tail;
406 if (ldu.d[1] >= 0.0) // round to integral
407 result = (ldu.d[1] + kTwoTo52) - kTwoTo52;
408 else
409 result = (ldu.d[1] - kTwoTo52) + kTwoTo52;
410
411 // make result "canonical" in to-nearest rounding (preserves value)
412 FESETENVD(0.0);
413
414 xd = ldu.d[0] + result;
415 ldu.d[1] = (ldu.d[0] - xd) + result;
416 ldu.d[0] = xd;
417
418 FESETENVD(fpenv); // restore caller's rounding
419 return (ldu.f); // deliver result
420}
421
422
423long double rintl( long double x )
424{
425 DblDbl ldu;
426 double fpenv;
427 Double fenv;
428 double result,xh,xt;
429 uint32_t rnddir;
430
431 FEGETENVD(fpenv);
432 FESETENVD(0.0);
433 fenv.f = fpenv;
434 rnddir = fenv.i[1] & FE_ALL_RND; // isolate round mode
435
436 if (rnddir == FE_TONEAREST) { // default rounding
437 ldu.f = x;
438
439 if ((ldu.i[0] & kExpMask) == kExpMask) {
440 FESETENVD(fpenv);
441 return (x);
442 }
443
444 if (fabs(ldu.d[1]) >= kTwoTo52) {
445 FESETENVD(fpenv);
446 return (x);
447 }
448
449 // binary point is within or to the left of x's bits; assume x is in
450 // canonical form for default rounding
451 xh = ldu.d[0]; // put x in canonical form
452 xt = ldu.d[1];
453
454 if (fabs(xh) < kTwoTo52) { // binary point in head
455 ldu.d[1] = 0.0; // result tail is zero
456 if (xh >= 0.0) // Rint(head)
457 result = (xh + kTwoTo52) - kTwoTo52;
458 else
459 result = (xh - kTwoTo52) + kTwoTo52;
460
461 // Fix up only false halfway cases
462 if ((fabs(result - xh) == 0.5) && (xt != 0.0))
463 result = (xt > 0.0) ? (xh + 0.5) : (xh - 0.5);
464
465 if (result == 0.0) { // preserve sign of zero
466 if (ldu.i[0] < kSignMask) {
467 result = 0.0;
468 }
469 else {
470 result = -0.0;
471 ldu.d[1] = result;
472 }
473 }
474 ldu.d[0] = result;
475 FESETENVD(fpenv);
476 return (ldu.f); // deliver result
477 } // [binary point in head]
478
479 // binary point in tail or between head and tail
480 if (xt >= 0.0) // round to integral
481 result = (xt + kTwoTo52) - kTwoTo52;
482 else
483 result = (xt - kTwoTo52) + kTwoTo52;
484
485 ldu.d[0] = xh + result; // make canonical
486 ldu.d[1] = xh - ldu.d[0] + result;
487 FESETENVD(fpenv);
488 return (ldu.f); // deliver result
489 } // [default rounding]
490
491 // non-default rounding
492 FESETENVD(fpenv);
493
494 if (rnddir == FE_TOWARDZERO) // round toward zero
495 return (truncl(x));
496 else if (rnddir == FE_UPWARD) // round upward
497 return (ceill(x));
498 else /* rnddir == FE_DOWNWARD */ // round downward
499 return (floorl(x));
500}
501
502long long int llrintl ( long double x )
503{
504 long double t;
505 long long int result;
506 double fenv;
507
508 if (unlikely(x != x))
509 {
510 feraiseexcept(FE_INVALID);
511 return LONG_LONG_MIN;
512 }
513
514 FEGETENVD(fenv);
515 t = rintl ( x );
516 FESETENVD(fenv);
517
518 if ( t < (long double)LONG_LONG_MIN )
519 {
520 feraiseexcept(FE_INVALID);
521 result = LONG_LONG_MIN;
522 }
523 else if ( t > (long double)LONG_LONG_MAX )
524 {
525 feraiseexcept(FE_INVALID);
526 result = LONG_LONG_MAX;
527 }
528 else if (t != x)
529 {
530 feraiseexcept(FE_INEXACT);
531 result = (long long int) t;
532 }
533 else
534 {
535 result = (long long int) t;
536 }
537
538 return result;
539}
540
541long double nearbyintl( long double x )
542{
543 return (rintl(x));
544}
545
546
547long double roundl( long double x )
548{
549 DblDbl ldu;
550 double xh, xt, xd, result, fpenv;
551
552 FEGETENVD(fpenv);
553 FESETENVD(0.0);
554
555 ldu.f = x;
556 xh = ldu.d[0];
557 xt = ldu.d[1];
558
559 if ((ldu.i[0] & kExpMask) == kExpMask) { // NaN or INF
560 FESETENVD(fpenv);
561 return (x);
562 }
563
564 if (fabs(xt) >= kTwoTo52) { // large integral x
565 FESETENVD(fpenv);
566 return (x);
567 }
568
569 // binary point is within or to the left of x's bits
570 FESETENVD(kTZ.f); // round toward zero
571
572 if (fabs(xh) < kTwoTo52) { // binary point in head
573 ldu.d[1] = 0.0; // result tail is zero
574 if (xh >= 0.0) {
575 xd = (xt + 0.5) + xh;
576 result = (xd + kTwoTo52) - kTwoTo52;
577 }
578 else {
579 xd = (xt - 0.5) + xh;
580 result = (xd - kTwoTo52) + kTwoTo52;
581 }
582
583 if (result == 0.0) { // preserve sign of zero
584 if (ldu.i[0] < kSignMask)
585 result = 0.0;
586 else {
587 result = -0.0;
588 ldu.d[1] = result;
589 }
590 }
591
592 ldu.d[0] = result;
593 FESETENVD(fpenv); // restore caller's env
594 return (ldu.f); // deliver result
595 } // [binary point in head]
596
597 // binary point in tail or between head and tail
598 if (xh > 0.0) {
599 xd = xt + 0.5;
600 FESETENVD(kDN.f);
601 }
602 else {
603 xd = xt - 0.5;
604 FESETENVD(kUP.f);
605 }
606
607 if (xd >= 0.0)
608 result = (xd + kTwoTo52) - kTwoTo52;
609 else
610 result = (xd - kTwoTo52) + kTwoTo52;
611
612 // make result "canonical" in to-nearest rounding (preserves value)
613 FESETENVD(0.0);
614 xd = ldu.d[0] + result;
615 ldu.d[1] = (ldu.d[0] - xd) + result;
616 ldu.d[0] = xd;
617 FESETENVD(fpenv);
618 return (ldu.f); // deliver result
619}
620
621long long int llroundl ( long double x )
622{
623 long double t;
624 long long int result;
625 double fenv;
626
627 if (unlikely(x != x))
628 {
629 feraiseexcept(FE_INVALID);
630 return LONG_LONG_MAX;
631 }
632
633 FEGETENVD(fenv);
634 t = roundl ( x );
635 FESETENVD(fenv);
636
637 if ( t < (long double)LONG_LONG_MIN )
638 {
639 feraiseexcept(FE_INVALID);
640 result = LONG_LONG_MIN;
641 }
642 else if ( t > (long double)LONG_LONG_MAX )
643 {
644 feraiseexcept(FE_INVALID);
645 result = LONG_LONG_MAX;
646 }
647 else if (t != x)
648 {
649 feraiseexcept(FE_INEXACT);
650 result = (long long int) t;
651 }
652 else
653 {
654 result = (long long int) t;
655 }
656
657 return result;
658}
659
660long double fdiml(long double x, long double y)
661{
662 double fenv;
663 long double result;
664
665 FEGETENVD(fenv);
666 FESETENVD(0.0);
667
668 if (unlikely((x != x) || (y != y)))
669 {
670 FESETENVD(fenv);
671 return ( x + y );
672 }
673 else if (x > y)
674 result = (x - y);
675 else
676 result = 0.0L;
677
678 FESETENVD(fenv);
679 return result;
680}
681
682long double fmaxl( long double x, long double y )
683{
684 DblDbl xDD,yDD;
685 double fpenv;
686
687 FEGETENVD(fpenv);
688 FESETENVD(0.0);
689
690 xDD.f = x;
691 yDD.f = y;
692
693 // filter out NaN arguments
694 if (yDD.d[0] != yDD.d[0]) {
695 FESETENVD(fpenv);
696 return (x);
697 }
698
699 if (xDD.d[0] != xDD.d[0]) {
700 FESETENVD(fpenv);
701 return (y);
702 }
703
704 if (y > x) {
705 FESETENVD(fpenv);
706 return (y);
707 }
708 else {
709 FESETENVD(fpenv);
710 return (x);
711 }
712}
713
714
715long double fminl( long double x, long double y )
716{
717 DblDbl xDD,yDD;
718 double fpenv;
719
720 FEGETENVD(fpenv);
721 FESETENVD(0.0);
722
723 xDD.f = x;
724 yDD.f = y;
725
726 // filter out NaN arguments
727 if (yDD.d[0] != yDD.d[0]) {
728 FESETENVD(fpenv);
729 return (x);
730 }
731
732 if (xDD.d[0] != xDD.d[0]) {
733 FESETENVD(fpenv);
734 return (y);
735 }
736
737 if (y < x) {
738 FESETENVD(fpenv);
739 return (y);
740 }
741 else {
742 FESETENVD(fpenv);
743 return (x);
744 }
745}
746
747long double fmal( long double x, long double y, long double z )
748{
749 DblDbl xDD, yDD, zDD;
750 double th, tm, tl, u;
751
752 xDD.f = x;
753 yDD.f = y;
754 zDD.f = z;
755
756 _d3mul( xDD.d[0], xDD.d[1], 0.0, yDD.d[0], yDD.d[1], 0.0, &th, &tm, &tl );
757 _d3add( th, tm, tl, zDD.d[0], zDD.d[1], 0.0, &xDD.d[0], &xDD.d[1], &u );
758
759 return xDD.f;
760}
761
762long double fabsl( long double x )
763{
764 DblDbl xDD;
765
766 xDD.f = x;
767
768 if (xDD.i[0] < kSignMask) // positive x
769 return (x);
770
771 else { // negative x
772 xDD.i[0] ^= kSignMask;
773 xDD.i[2] ^= kSignMask;
774 return (xDD.f);
775 }
776}
777
778
779long double copysignl( long double x, long double y )
780{
781 DblDbl xDD,yDD;
782
783 xDD.f = x;
784 yDD.f = y;
785
786 if ((xDD.i[0] & kSignMask) == (yDD.i[0] & kSignMask))
787 return (x);
788
789 else { // signs different
790 xDD.i[0] ^= kSignMask;
791 xDD.i[2] ^= kSignMask;
792 return (xDD.f);
793 }
794}
795
796
797long rinttoll( long double x )
798{
799 DblDbl ldu;
800 uint32_t expx, lowbit;
801 long temp;
802 double fpenv;
803
804 ldu.f = x;
805 expx = (ldu.i[0] >> 20) & 0x7ffu; // biased exp
806 lowbit = ldu.i[1] & 1u; // low head sig bit
807
808 // if x is roughly in range of long format and tail component is nonzero,
809 // adjust its head component to contain proper value and sticky information
810 // for conversion to long.
811
812 if ((expx < 0x41fu) && (expx > 0x3fdu)) { // 2^32 > |x| > (0.5 - 1 ulp)
813 if ((lowbit == 0) && (ldu.d[1] != 0.0)) {
814 uint32_t signh = (ldu.i[0] >> 31) & 1u; // sign of head
815 uint32_t signt = (ldu.i[2] >> 31) & 1u; // sign of tail
816 if (signh == signt) // same sign
817 ldu.i[1] |= 1u;
818 else {
819 ldu.i[1] -= 1u; // unlike signs
820 if (ldu.i[1] == 0xffffffffu)
821 ldu.i[0] -= 1u;
822 }
823 }
824 }
825 FEGETENVD(fpenv);
826 FESETENVD(0.0);
827 temp = rinttol(ldu.d[0]);
828 FESETENVD(fpenv);
829
830// return (rinttol(ldu.d[0]));
831 return (temp);
832}
833
834// These work just as well for the LP64 ABI
835long int lroundl( long double x )
836{
837 long double t;
838 long int result;
839 double fenv;
840
841 if (unlikely(x != x))
842 {
843 feraiseexcept(FE_INVALID);
844 return LONG_MAX;
845 }
846
847 FEGETENVD(fenv);
848 t = roundl ( x );
849 FESETENVD(fenv);
850
851 if ( t < (long double)LONG_MIN )
852 {
853 feraiseexcept(FE_INVALID);
854 result = LONG_MIN;
855 }
856 else if ( t > (long double)LONG_MAX )
857 {
858 feraiseexcept(FE_INVALID);
859 result = LONG_MAX;
860 }
861 else if (t != x)
862 {
863 feraiseexcept(FE_INEXACT);
864 result = (long int) t;
865 }
866 else
867 {
868 result = (long int) t;
869 }
870
871 return result;
872}
873
874long double ldexpl(long double x, int n)
875{
876 int m;
877
878 // Clip n
879 if (unlikely(n > 2097))
880 m = 2098;
881 else if (unlikely(n < -2098))
882 m = -2099;
883 else
884 m = n;
885
886 return scalbnl(x, m);
887}
888
889long double frexpl(long double x, int *exponent)
890{
891 DblDbl ldu;
892 double xHead;
893 double fenv;
894 long double result;
895
896 ldu.f = x;
897 xHead = ldu.d[0];
898
899 FEGETENVD(fenv);
900 FESETENVD(0.0);
901
902 switch (__fpclassifyd(xHead))
903 {
904 //case FP_SNAN:
905 //case FP_QNAN:
906 case FP_NAN:
907 case FP_INFINITE:
908 FESETENVD(fenv);
909 return x;
910 default:
911 *exponent = (x == 0.0) ? 0 : (int)(1.0 + logbl(x));
912 result = scalbnl(x, - (*exponent));
913 FESETENVD(fenv);
914 return result;
915 }
916}
917//
918//
919//relop relationl(long double x, long double y)
920//{
921// double fenv;
922// relop result;
923//
924// fenv = __setflm(0.0);
925//
926// if (( x != x ) || ( y != y ))
927// result = UNORDERED;
928// else if (x == y)
929// result = EQUALTO;
930// else if (x < y)
931// result = LESSTHAN;
932// else
933// result = GREATERTHAN;
934//
935// __setflm(fenv);
936// return result;
937//}
938
939long double fmodl(long double xx, long double yy)
940{
941 int iclx,icly; /* classify results of x,y */
942 int32_t iscx,iscy,idiff; /* logb values and difference */
943 int i; /* loop variable */
944 long double absy,x1,y1,z; /* local floating-point variables */
945 long double rslt;
946 fenv_t OldEnv;
947 hexdouble OldEnvironment;
948 int newexc;
949 DblDbl xDD, yDD;
950
951 xDD.f = xx;
952 yDD.f = yy;
953
954 FEGETENVD( OldEnvironment.d );
955 FESETENVD( 0.0 );
956 OldEnv = OldEnvironment.i.lo;
957
958 iclx = __fpclassifyd(xDD.d[0]);
959 icly = __fpclassifyd(yDD.d[0]);
960 if (likely((iclx & icly) >= FP_NORMAL)) { /* x,y both nonzero finite case */
961 x1 = fabsl(xx); /* work with absolute values */
962 absy = fabsl(yy);
963 if (absy > x1) {
964 rslt = xx; /* trivial case */
965 goto ret;
966 }
967 else { /* nontrivial case requires reduction */
968 iscx = (int32_t) logbl(x1); /* get binary exponents of |x| and |y| */
969 iscy = (int32_t) logbl(absy);
970 idiff = iscx - iscy; /* exponent difference */
971 if (idiff != 0) { /* exponent of x1 > exponent of y1 */
972 y1 = scalbnl(absy,-iscy); /* scale |y| to unit binade */
973 x1 = scalbnl(x1,-iscx); /* ditto for |x| */
974 for (i = idiff; i != 0; i--) { /* begin remainder loop */
975 if ((z = x1 - y1) >= 0) { /* nonzero remainder step result */
976 x1 = z; /* update remainder (x1) */
977 }
978 x1 += x1; /* shift (by doubling) remainder */
979 } /* end of remainder loop */
980 x1 = scalbnl(x1,iscy); /* scale result to binade of |y| */
981 } /* remainder exponent >= exponent of y */
982 if (x1 >= absy) { /* last step to obtain modulus */
983 x1 -= absy;
984 }
985 } /* x1 is |result| */
986 if (xx < 0.0)
987 x1 = -x1; /* modulus if x is negative */
988 rslt = x1;
989 goto ret;
990 } /* end of x,y both nonzero finite case */
991 else if ((iclx <= FP_QNAN) || (icly <= FP_QNAN)) {
992 rslt = xx+yy; /* at least one NaN operand */
993 goto ret;
994 }
995 else if ((iclx == FP_INFINITE)||(icly == FP_ZERO)) { /* invalid result */
996 rslt = nanl(REM_NAN);
997 OldEnvironment.i.lo |= SET_INVALID;
998 FESETENVD ( OldEnvironment.d );
999 goto ret;
1000 }
1001 else /* trivial cases (finite MOD infinite */
1002 rslt = xx; /* or zero REM nonzero) with *quo = 0 */
1003 ret:
1004 FEGETENVD (OldEnvironment.d );
1005 newexc = OldEnvironment.i.lo & FE_ALL_EXCEPT;
1006 OldEnvironment.i.lo = OldEnv;
1007 if ((newexc & FE_INVALID) != 0)
1008 OldEnvironment.i.lo |= SET_INVALID;
1009 OldEnvironment.i.lo |= newexc & ( FE_INEXACT | FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW );
1010 FESETENVD (OldEnvironment.d );
1011 return rslt;
1012}
1013
1014#warning remquol: cannot gaurantee exact result!
1015static const hexdouble Huge = HEXDOUBLE(0x7ff00000, 0x00000000);
1016static const hexdouble HugeHalved = HEXDOUBLE(0x7fe00000, 0x00000000);
1017long double remquol ( long double x, long double y, int *quo)
1018{
1019 int iclx,icly; /* classify results of x,y */
1020 int32_t iquo; /* low 32 bits of integral quotient */
1021 int32_t iscx, iscy, idiff; /* logb values and difference */
1022 int i; /* loop variable */
1023 long double absy,x1,y1,z; /* local floating-point variables */
1024 long double rslt;
1025 fenv_t OldEnv;
1026 hexdouble OldEnvironment;
1027 int newexc;
1028
1029 FEGETENVD ( OldEnvironment.d );
1030 FESETENVD ( 0.0 );
1031 __NOOP;
1032 __NOOP;
1033
1034 OldEnv = OldEnvironment.i.lo;
1035
1036 *quo = 0; /* initialize quotient result */
1037 iclx = fpclassify(x);
1038 icly = fpclassify(y);
1039 if (likely((iclx & icly) >= FP_NORMAL)) { /* x,y both nonzero finite case */
1040 x1 = fabsl(x); /* work with absolute values */
1041 absy = fabsl(y);
1042 iquo = 0; /* zero local quotient */
1043 iscx = (int32_t) logbl(x1); /* get binary exponents */
1044 iscy = (int32_t) logbl(absy);
1045 idiff = iscx - iscy; /* exponent difference */
1046 if (idiff >= 0) { /* exponent of x1 >= exponent of y1 */
1047 if (idiff != 0) { /* exponent of x1 > exponent of y1 */
1048 y1 = scalbnl(absy,-iscy); /* scale |y| to unit binade */
1049 x1 = scalbnl(x1,-iscx); /* ditto for |x| */
1050 for (i = idiff; i != 0; i--) { /* begin remainder loop */
1051 if ((z = x1 - y1) >= 0) { /* nonzero remainder step result */
1052 x1 = z; /* update remainder (x1) */
1053 iquo += 1; /* increment quotient */
1054 }
1055 iquo += iquo; /* shift quotient left one bit */
1056 x1 += x1; /* shift (double) remainder */
1057 } /* end of remainder loop */
1058 x1 = scalbnl(x1,iscy); /* scale remainder to binade of |y| */
1059 } /* remainder has exponent <= exponent of y */
1060 if (x1 >= absy) { /* last remainder step */
1061 x1 -= absy;
1062 iquo +=1;
1063 } /* end of last remainder step */
1064 } /* remainder (x1) has smaller exponent than y */
1065 if (likely( x1 < HugeHalved.d ))
1066 z = scalbnl( x1, 1 ); // z = x1 + x1; /* double remainder, without overflow */
1067 else
1068 z = Huge.d;
1069 if ((z > absy) || ((z == absy) && ((iquo & 1) != 0))) {
1070 x1 -= absy; /* final remainder correction */
1071 iquo += 1;
1072 }
1073 if (x < 0.0)
1074 x1 = -x1; /* remainder if x is negative */
1075 iquo &= 0x0000007f; /* retain low 7 bits of integer quotient */
1076 if ((signbit(x) ^ signbit(y)) != 0) /* take care of sign of quotient */
1077 iquo = -iquo;
1078 *quo = iquo; /* deliver quotient result */
1079 rslt = x1;
1080 goto ret;
1081 } /* end of x,y both nonzero finite case */
1082 else if ((iclx <= FP_QNAN) || (icly <= FP_QNAN)) {
1083 rslt = x+y; /* at least one NaN operand */
1084 goto ret;
1085 }
1086 else if ((iclx == FP_INFINITE)||(icly == FP_ZERO)) { /* invalid result */
1087 rslt = nan("9");
1088 OldEnvironment.i.lo |= SET_INVALID;
1089 FESETENVD_GRP( OldEnvironment.d );
1090 goto ret;
1091 }
1092 else /* trivial cases (finite REM infinite */
1093 rslt = x; /* or zero REM nonzero) with *quo = 0 */
1094 ret:
1095 FEGETENVD_GRP( OldEnvironment.d );
1096 newexc = OldEnvironment.i.lo & FE_ALL_EXCEPT;
1097 OldEnvironment.i.lo = OldEnv;
1098 if ((newexc & FE_INVALID) != 0)
1099 OldEnvironment.i.lo |= SET_INVALID;
1100 OldEnvironment.i.lo |= newexc & ( FE_INEXACT | FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW );
1101 FESETENVD_GRP( OldEnvironment.d );
1102 return rslt;
1103}
1104
1105long double remainderl(long double x, long double y)
1106{
1107 int quo;
1108 return (remquol(x, y, &quo));
1109}
1110
1111long double modfl ( long double x, long double *iptr )
1112{
1113 DblDbl ldu;
1114 double xh, xt, frach, fract, inth, intt, fpenv;
1115 long double l;
1116
1117 FEGETENVD(fpenv);
1118 FESETENVD(0.0);
1119
1120 ldu.f = x;
1121 xh = ldu.d[0];
1122 xt = ldu.d[1];
1123
1124 frach = modf( xh, &inth );
1125 fract = modf( xt, &intt );
1126
1127 *iptr = (long double)(inth + intt);
1128 l = ((long double)frach) + ((long double)fract);
1129 if (l >= 1.0L)
1130 {
1131 l = l - 1.0L;
1132 *iptr = *iptr + 1.0L;
1133 }
1134
1135
1136 FESETENVD(fpenv);
1137 return l;
1138}
1139
1140long double nextafterl(long double x, long double y)
1141{
1142 static const hexdouble EPSILON = HEXDOUBLE(0x00000000, 0x00000001);
1143 static const hexdouble PosBigHead = HEXDOUBLE(0x7fefffff, 0xffffffff);
1144 static const hexdouble PosBigTail = HEXDOUBLE(0x7c9fffff, 0xffffffff);
1145
1146 DblDbl ldu;
1147 double fpenv;
1148
1149 if (unlikely( ( x != x ) || ( y != y ) )) // one of the arguments is a NaN
1150 return x + y;
1151 else if ( x == y )
1152 {
1153 if ( ( x == 0.0L ) && ( y == 0.0L )) // (*0, -0)-> -0, (*0, +0)-> +0 i.e. "y"
1154 return y;
1155 else
1156 return x;
1157 }
1158 else if (unlikely( isinf( x ) )) // N.B. y is not equal to (infinite) x hence is finite
1159 {
1160 ldu.d[0] = PosBigHead.d;
1161 ldu.d[1] = PosBigTail.d;
1162
1163 return copysignl( ldu.f, x );
1164 }
1165 else if (x == 0.0L)
1166 {
1167 ldu.d[0] = EPSILON.d;
1168 ldu.d[1] = 0.0L;
1169
1170 return copysignl( ldu.f, y );
1171 }
1172 else if ( ( ( x < 0.0 ) && ( x < y ) ) || ( ( x > 0.0 ) && ( x > y ) ) )
1173 {
1174 FEGETENVD(fpenv);
1175 FESETENVD(kTZ.f);
1176 x = (x < 0.0L) ? x + (long double)EPSILON.d : x - (long double)EPSILON.d;
1177 FESETENVD(fpenv);
1178 return x;
1179 }
1180 else if ( ( x < 0.0 ) && ( x > y ) )
1181 {
1182 FEGETENVD(fpenv);
1183 FESETENVD(kDN.f);
1184 x = x - (long double)EPSILON.d;
1185 FESETENVD(fpenv);
1186 return x;
1187 }
1188 else // ( ( x > 0.0 ) && ( x < y ) )
1189 {
1190 FEGETENVD(fpenv);
1191 FESETENVD(kUP.f);
1192 x = x + (long double)EPSILON.d;
1193 FESETENVD(fpenv);
1194 return x;
1195 }
1196}
1197
1198float nexttowardf(float x, long double y)
1199{
1200 DblDbl ldu;
1201 double yh, yt, yd, fpenv;
1202
1203 if ((long double)x == y)
1204 return (float)y; // 7.12.11.4 requires this behavior
1205
1206 FEGETENVD(fpenv);
1207 FESETENVD(0.0);
1208
1209 ldu.f = y;
1210 yh = ldu.d[0];
1211 yt = ldu.d[1];
1212
1213 yd = yh + yt;
1214 FESETENVD(fpenv);
1215
1216 return nextafterf(x, (float)yd);
1217}
1218
1219double nexttoward(double x, long double y)
1220{
1221 DblDbl ldu;
1222 double yh, yt, yd, fpenv;
1223
1224 if ((long double)x == y)
1225 return (double)y; // 7.12.11.4 requires this behavior
1226
1227 FEGETENVD(fpenv);
1228 FESETENVD(0.0);
1229
1230 ldu.f = y;
1231 yh = ldu.d[0];
1232 yt = ldu.d[1];
1233
1234 yd = yh + yt;
1235 FESETENVD(fpenv);
1236
1237 return nextafter(x, y);
1238}
1239
1240long double nexttowardl(long double x, long double y)
1241{
1242 if (x == y)
1243 return y; // 7.12.11.4 requires this behavior
1244 else
1245 return nextafterl( x, y );
1246}
1247#endif
1248