home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2337 / nextafter.c next >
Encoding:
C/C++ Source or Header  |  1990-12-28  |  1.7 KB  |  99 lines
  1. /*
  2. ** This file is part of the alternative 80386 math library and is
  3. ** covered by the GNU General Public license with my modification
  4. ** as noted in the README file that accompanied this file.
  5. **
  6. ** Copyright 1990 G. Geers
  7. **
  8. ** A mix of C and assembler - well I've got the functions so I might 
  9. ** as well use them!
  10. **
  11. */
  12.  
  13. #include "fpumath.h"
  14.  
  15. asm(".align 4");
  16. asm(".Lulp:");
  17. asm(".double 1.1102230246251565e-16");
  18.  
  19. asm(".align 4");
  20. asm(".Lulpup:");
  21. asm(".double 2.2204460492503131e-16");
  22.  
  23. double
  24. nextafter(x, y)
  25. double x, y;
  26. {
  27.     asm("subl $8, %esp");
  28.  
  29.     if (isnan(x) || isnan(y))
  30.         return(quiet_nan(1.0));
  31.  
  32.     if (isinf(x))
  33.         if (y > x)
  34.             return(-max_normal());
  35.         else
  36.             if (y < x)
  37.                 return(max_normal());
  38.  
  39.     if (x == 0.0)
  40.         if (y > 0.0)
  41.             return(min_subnormal());
  42.         else
  43.             return(-min_subnormal());
  44.  
  45.     if (isnormal(x)) {
  46.         if ((x == min_normal()) && y < x)
  47.             return(max_subnormal());
  48.  
  49.         if ((x == max_normal()) && y > x)
  50.             return(infinity());
  51.  
  52.         if ((x == -max_normal()) && y < x)
  53.             return(-infinity());
  54.  
  55.         asm("movl 12(%ebp), %eax");
  56.         asm("andl $0x7ff00000, %eax");
  57.         asm("movl %eax, -12(%ebp)");
  58.         asm("movl $0x0, -16(%ebp)");
  59.         asm("fincstp");
  60.  
  61.         if (fabs(x) <= 2.0 && y < x)
  62.             asm("fldl .Lulp");
  63.         else
  64.             asm("fldl .Lulpup");
  65.  
  66.         asm("fldl -16(%ebp)");
  67.         asm("fmulp");
  68.  
  69.         if (y > x) {
  70.             asm("fldl 8(%ebp)");
  71.             asm("faddp");
  72.             asm("leave");
  73.             asm("ret");
  74.         }
  75.         if (y < x) {
  76.             asm("fldl 8(%ebp)");
  77.             asm("fsubp");
  78.             asm("leave");
  79.             asm("ret");
  80.         }
  81.     }
  82.     else
  83.     if (issubnormal(x)) {
  84.         if ((x == max_subnormal()) && y > x)
  85.             return(min_normal());
  86.  
  87.         if ((x == -max_subnormal()) && y < x)
  88.             return(-min_normal());
  89.  
  90.         if (y > x) 
  91.             return(x + min_subnormal());
  92.  
  93.         if (y < x)
  94.             return(x - min_subnormal());
  95.     }
  96.     
  97.     return(x);
  98. }
  99.