This is mpmath.c in view mode; [Download] [Up]
/****************************************************************
Copyright 1990 - 1995 by AT&T Bell Laboratories.
Permission to use, copy, modify, and distribute this software
and its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the names of AT&T Bell Laboratories or
any of its entities not be used in advertising or publicity
pertaining to distribution of the software without specific,
written prior permission.
AT&T disclaims all warranties with regard to this software,
including all implied warranties of merchantability and fitness.
In no event shall AT&T be liable for any special, indirect or
consequential damages or any damages whatsoever resulting from
loss of use, data or profits, whether in an action of contract,
negligence or other tortious action, arising out of or in
connection with the use or performance of this software.
****************************************************************/
#define EXTERN extern
#include "mpd.h" /* Includes "site.h" via "mp.h" */
/**********************************************************
The following is by John Hobby
**********************************************************/
#ifndef FIXPT
/* These replacements for takefraction, makefraction, takescaled, makescaled
run about 3 to 11 times faster than the standard versions on modern machines
that have fast hardware for double-precision floating point. They should
produce approximately correct results on all machines and agree exactly
with the standard versions on machines that satisfy the following conditions:
1. Doubles must have at least 46 mantissa bits; i.e., numbers expressible
as n*2^k with abs(n)<2^46 should be representable.
2. The following should hold for addition, subtraction, and multiplcation but
not necessarily for division:
A. If the true answer is between two representable numbers, the computed
answer must be one of them.
B. When the true answer is representable, this must be the computed result.
3. Dividing one double by another should always produce a relative error of
at most one part in 2^46. (This is why the mantissa requirement is
46 bits instead of 45 bits.)
3. In the absence of overflow, double-to-integer conversion should truncate
toward zero and do this in an exact fashion.
4. Integer-to-double convesion should produce exact results.
5. Dividing one power of two by another should yield an exact result.
6. ASCII to double conversion should be exact for integer values.
7. Integer arithmetic must be done in the two's-complement system.
*/
#define ELGORDO 0x7fffffff
#define TWEXP31 2147483648.0
#define TWEXP28 268435456.0
#define TWEXP16 65536.0
#define TWEXP_16 (1.0/65536.0)
#define TWEXP_28 (1.0/268435456.0)
integer ztakefraction(p,q) /* Approximate p*q/2^28 */
integer p,q;
{ register double d;
register integer i;
d = (double)p * (double)q * TWEXP_28;
if ((p^q) >= 0) {
d += 0.5;
if (d>=TWEXP31) {
if (d!=TWEXP31 || (((p&077777)*(q&077777))&040000)==0)
aritherror = true;
return ELGORDO;
}
i = (integer) d;
if (d==i && (((p&077777)*(q&077777))&040000)!=0) --i;
} else {
d -= 0.5;
if (d<= -TWEXP31) {
if (d!= -TWEXP31 || ((-(p&077777)*(q&077777))&040000)==0)
aritherror = true;
return -ELGORDO;
}
i = (integer) d;
if (d==i && ((-(p&077777)*(q&077777))&040000)!=0) ++i;
}
return i;
}
integer ztakescaled(p,q) /* Approximate p*q/2^16 */
integer p,q;
{ register double d;
register integer i;
d = (double)p * (double)q * TWEXP_16;
if ((p^q) >= 0) {
d += 0.5;
if (d>=TWEXP31) {
if (d!=TWEXP31 || (((p&077777)*(q&077777))&040000)==0)
aritherror = true;
return ELGORDO;
}
i = (integer) d;
if (d==i && (((p&077777)*(q&077777))&040000)!=0) --i;
} else {
d -= 0.5;
if (d<= -TWEXP31) {
if (d!= -TWEXP31 || ((-(p&077777)*(q&077777))&040000)==0)
aritherror = true;
return -ELGORDO;
}
i = (integer) d;
if (d==i && ((-(p&077777)*(q&077777))&040000)!=0) ++i;
}
return i;
}
/* Note that d cannot exactly equal TWEXP31 when the overflow test is made
because the exact value of p/q cannot be strictly between (2^31-1)/2^28
and 8/1. No pair of integers less than 2^31 has such a ratio.
*/
integer zmakefraction(p,q) /* Approximate 2^28*p/q */
integer p,q;
{ register double d;
register integer i;
#ifdef DEBUG
if (q==0) confusion(47);
#endif /* DEBUG */
d = TWEXP28 * (double)p /(double)q;
if ((p^q) >= 0) {
d += 0.5;
if (d>=TWEXP31) {aritherror=true; return ELGORDO;}
i = (integer) d;
if (d==i && ( ((q>0 ? -q : q)&077777)
* (((i&037777)<<1)-1) & 04000)!=0) --i;
} else {
d -= 0.5;
if (d<= -TWEXP31) {aritherror=true; return -ELGORDO;}
i = (integer) d;
if (d==i && ( ((q>0 ? q : -q)&077777)
* (((i&037777)<<1)+1) & 04000)!=0) ++i;
}
return i;
}
/* Note that d cannot exactly equal TWEXP31 when the overflow test is made
because the exact value of p/q cannot be strictly between (2^31-1)/2^16
and 2^15/1. No pair of integers less than 2^31 has such a ratio.
*/
integer zmakescaled(p,q) /* Approximate 2^16*p/q */
integer p,q;
{ register double d;
register integer i;
#ifdef DEBUG
if (q==0) confusion(47);
#endif /* DEBUG */
d = TWEXP16 * (double)p /(double)q;
if ((p^q) >= 0) {
d += 0.5;
if (d>=TWEXP31) {aritherror=true; return ELGORDO;}
i = (integer) d;
if (d==i && ( ((q>0 ? -q : q)&077777)
* (((i&037777)<<1)-1) & 04000)!=0) --i;
} else {
d -= 0.5;
if (d<= -TWEXP31) {aritherror=true; return -ELGORDO;}
i = (integer) d;
if (d==i && ( ((q>0 ? q : -q)&077777)
* (((i&037777)<<1)+1) & 04000)!=0) ++i;
}
return i;
}
#endif /* not FIXPT */
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.