home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / mitsch75.zip / scheme-7_5_17-src.zip / scheme-7.5.17 / src / microcode / missing.c < prev    next >
C/C++ Source or Header  |  2000-12-05  |  5KB  |  257 lines
  1. /* -*-C-*-
  2.  
  3. $Id: missing.c,v 9.33 2000/12/05 21:23:45 cph Exp $
  4.  
  5. Copyright (c) 1987-2000 Massachusetts Institute of Technology
  6.  
  7. This program is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU General Public License as published by
  9. the Free Software Foundation; either version 2 of the License, or (at
  10. your option) any later version.
  11.  
  12. This program is distributed in the hope that it will be useful, but
  13. WITHOUT ANY WARRANTY; without even the implied warranty of
  14. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15. General Public License for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with this program; if not, write to the Free Software
  19. Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20. */
  21.  
  22. /* This file contains utilities potentially missing from the math library. */
  23.  
  24. #include "config.h"
  25.  
  26. #ifndef HAVE_FREXP
  27.  
  28. double
  29. DEFUN (frexp, (value, eptr),
  30.        double value
  31.        AND int * eptr)
  32. {
  33.   register double x = ((value < 0) ? (-value) : value);
  34.   int e = 0;
  35.   if (x >= 1)
  36.     {
  37.       while (1)
  38.     {
  39.       if (x > 2)
  40.         {
  41.           register double xr = (x / 2);
  42.           register double r = 2;
  43.           register int n = 1;
  44.           while (xr >= r)
  45.         {
  46.           /* ((xr == (x / r)) && (xr >= r) && (x >= (r * r))) */
  47.           xr /= r;
  48.           /* ((xr == (x / (r * r))) && (xr >= 1)) */
  49.           r *= r;
  50.           /* ((xr == (x / r)) && (x >= r)) */
  51.           n += n;
  52.         }
  53.           /* ((xr >= 1) && (xr < r)) */
  54.           x = xr;
  55.           e += n;
  56.         }
  57.       else if (x < 2)
  58.         {
  59.           x /= 2;
  60.           e += 1;
  61.           break;
  62.         }
  63.       else
  64.         {
  65.           x /= 4;
  66.           e += 2;
  67.           break;
  68.         }
  69.     }
  70.     }
  71.   else if ((x > 0) && (x < 0.5))
  72.     {
  73.       while (1)
  74.     {
  75.       if (x < 0.25)
  76.         {
  77.           register double xr = (x * 2);
  78.           register double r = 0.5;
  79.           register int n = 1;
  80.           /* ((xr == (x / r)) && (xr < 0.5) && (x < (r / 2))) */
  81.           while (xr < (r / 2))
  82.         {
  83.           /* ((xr < (r / 2)) && (x < ((r * r) / 2))) */
  84.           xr /= r;
  85.           /* ((xr == (x / (r * r))) && (xr < 0.5)) */
  86.           r *= r;
  87.           /* ((xr == (x / r)) && (x < (r / 2))) */
  88.           n += n;
  89.         }
  90.           /* ((xr >= (r / 2)) && (xr < 0.5)) */
  91.           x = xr;
  92.           e -= n;
  93.         }
  94.       else
  95.         {
  96.           x *= 2;
  97.           e -= 1;
  98.           break;
  99.         }
  100.     }
  101.     }
  102.   (*eptr) = e;
  103.   return ((value < 0) ? (-x) : x);
  104. }
  105.  
  106. double
  107. DEFUN (ldexp, (value, exponent),
  108.        double value
  109.        AND int exponent)
  110. {
  111.   register double x = value;
  112.   register int e = exponent;
  113.   register double r = 2;
  114.   if (e > 0)
  115.     {
  116.       if (e == 1)
  117.     return (x * 2);
  118.       while (1)
  119.     {
  120.       if ((e % 2) != 0)
  121.         x *= r;
  122.       e /= 2;
  123.       if (e == 1)
  124.         return (x * r * r);
  125.       r *= r;
  126.     }
  127.     }
  128.   else if (e < 0)
  129.     {
  130.       e = (-e);
  131.       if (e == 1)
  132.     return (x / 2);
  133.       while (1)
  134.     {
  135.       if ((e % 2) != 0)
  136.         x /= r;
  137.       e /= 2;
  138.       if (e == 1)
  139.         return ((x / r) / r);
  140.       r *= r;
  141.     }
  142.     }
  143.   else
  144.     return (x);
  145. }
  146.  
  147. #endif /* not HAVE_FREXP */
  148.  
  149. #ifndef HAVE_MODF
  150.  
  151. double
  152. DEFUN (modf, (value, iptr),
  153.        double value
  154.        AND double * iptr)
  155. {
  156.   int exponent;
  157.   double significand = (frexp (value, (&exponent)));
  158.   if ((significand == 0) || (exponent <= 0))
  159.     {
  160.       (*iptr) = 0;
  161.       return (value);
  162.     }
  163.   {
  164.     register double s =
  165.       ((((significand < 0) ? (-significand) : significand) * 2) - 1);
  166.     register int e = (exponent - 1);
  167.     register double n = 1;
  168.     while (1)
  169.       {
  170.     if (e == 0)
  171.       break;
  172.     s *= 2;
  173.     e -= 1;
  174.     n *= 2;
  175.     if (s >= 1)
  176.       {
  177.         s -= 1;
  178.         n += 1;
  179.         if (s <= 0)
  180.           {
  181.         /* Multiply n by 2^e */
  182.         register double b = 2;
  183.         if (e == 0)
  184.           break;
  185.         while (1)
  186.           {
  187.             if ((e % 2) == 1)
  188.               {
  189.             n *= b;
  190.             if (e == 1)
  191.               break;
  192.             e -= 1;
  193.               }
  194.             b *= b;
  195.             e /= 2;
  196.           }
  197.         break;
  198.           }
  199.       }
  200.       }
  201.     if (significand < 0)
  202.     {
  203.       (*iptr) = (- n);
  204.       return (- s);
  205.     }
  206.     else
  207.     {
  208.       (*iptr) = n;
  209.       return (s);
  210.     }
  211.   }
  212. }
  213.  
  214. #endif /* not HAVE_MODF */
  215.  
  216. #ifndef HAVE_FLOOR
  217.  
  218. double
  219. DEFUN (floor, (x), double x)
  220. {
  221.   double iptr;
  222.   double fraction = (modf (x, (&iptr)));
  223.   return ((fraction < 0) ? (iptr - 1) : iptr);
  224. }
  225.  
  226. double
  227. DEFUN (ceil, (x), double x)
  228. {
  229.   double iptr;
  230.   double fraction = (modf (x, (&iptr)));
  231.   return ((fraction > 0) ? (iptr + 1) : iptr);
  232. }
  233.  
  234. #endif /* not HAVE_FLOOR */
  235.  
  236. #ifdef DEBUG_MISSING
  237.  
  238. #include <stdio.h>
  239.  
  240. main ()
  241. {
  242.   double input;
  243.   double output;
  244.   int exponent;
  245.   while (1)
  246.     {
  247.       printf ("Number -> ");
  248.       scanf ("%F", (&input));
  249.       output = (frexp (input, (&exponent)));
  250.       printf ("Input = %G; Output = %G; Exponent = %d\n",
  251.           input, output, exponent);
  252.       printf ("Result = %G\n", (ldexp (output, exponent)));
  253.     }
  254. }
  255.  
  256. #endif /* DEBUG_MISSING */
  257.