"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "scl.c" between
scm-5f2.zip and scm-5f3.zip

About: SCM is a Scheme Language Interpreter.

scl.c  (scm-5f2):scl.c  (scm-5f3)
skipping to change at line 43 skipping to change at line 43
# endif # endif
static sizet idbl2str P((double f, char *a)); static sizet idbl2str P((double f, char *a));
static sizet pdbl2str P((double f, char *a, sizet ch)); static sizet pdbl2str P((double f, char *a, sizet ch));
static sizet iflo2str P((SCM flt, char *str)); static sizet iflo2str P((SCM flt, char *str));
static void safe_add_1 P((double f, double *fsum)); static void safe_add_1 P((double f, double *fsum));
static long scm_twos_power P((SCM n)); static long scm_twos_power P((SCM n));
static double mantexp2dbl P((SCM manstr, int expo)); static double mantexp2dbl P((SCM manstr, int expo));
static double pmantexp2dbl P((SCM bmant, int expo)); static double pmantexp2dbl P((SCM bmant, int expo));
static SCM ex2in P((SCM z)); static SCM ex2in P((SCM z));
static SCM ilog P((unsigned long m, SCM b, SCM k, unsigned long *n)); static SCM ilog P((unsigned long m, SCM b, SCM k, unsigned long *n));
static double dpows5[23];
static char s_makrect[] = "make-rectangular", s_makpolar[] = "make-polar", static char s_makrect[] = "make-rectangular", s_makpolar[] = "make-polar",
s_magnitude[] = "magnitude", s_angle[] = "angle", s_magnitude[] = "magnitude", s_angle[] = "angle",
s_real_part[] = "real-part", s_imag_part[] = "imag-part", s_real_part[] = "real-part", s_imag_part[] = "imag-part",
s_in2ex[] = "inexact->exact",s_ex2in[] = "exact->inexact"; s_in2ex[] = "inexact->exact",s_ex2in[] = "exact->inexact";
static char s_expt[] = "real-expt", s_atan2[] = "$atan2"; static char s_expt[] = "real-expt", s_atan2[] = "$atan2";
#endif #endif
static char s_memv[] = "memv", s_assv[] = "assv"; static char s_memv[] = "memv", s_assv[] = "assv";
skipping to change at line 72 skipping to change at line 73
static char s_max[] = "max", s_min[] = "min"; static char s_max[] = "max", s_min[] = "min";
char s_sum[] = "+", s_difference[] = "-", s_product[] = "*", char s_sum[] = "+", s_difference[] = "-", s_product[] = "*",
s_divide[] = "/"; s_divide[] = "/";
static char s_number2string[] = "number->string", static char s_number2string[] = "number->string",
s_str2number[] = "string->number"; s_str2number[] = "string->number";
static char s_list_tail[] = "list-tail"; static char s_list_tail[] = "list-tail";
static char s_str2list[] = "string->list"; static char s_str2list[] = "string->list";
static char s_st_copy[] = "string-copy", s_st_fill[] = "string-fill!"; static char s_st_copy[] = "string-copy", s_st_fill[] = "string-fill!";
static char s_vect2list[] = "vector->list", s_ve_fill[] = "vector-fill!"; static char s_vect2list[] = "vector->list", s_ve_fill[] = "vector-fill!";
static char str_inf0[] = "inf.0", str_nan0[] = "nan.0", str_0[] = "0.0"; static char str_inf0[] = "inf.0", str_nan0[] = "nan.0", str_0[] = "0.";
static char s_intexpt[] = "integer-expt", s_cintlog[] = "ceiling-integer-log"; static char s_intexpt[] = "integer-expt", s_cintlog[] = "ceiling-integer-log";
static char s_dfloat_parts[] = "double-float-parts"; static char s_dfloat_parts[] = "double-float-parts";
#define s_intlog (&s_cintlog[8]) #define s_intlog (&s_cintlog[8])
/*** NUMBERS -> STRINGS ***/ /*** NUMBERS -> STRINGS ***/
#ifdef FLOATS #ifdef FLOATS
static int dbl_mant_dig = 0; static int dbl_mant_dig = 0;
static double max_dbl_int; /* Integers less than or equal to max_dbl_int static double max_dbl_int; /* Integers less than or equal to max_dbl_int
are representable exactly as doubles. */ are representable exactly as doubles. */
skipping to change at line 169 skipping to change at line 170
ndig -= DBL_MIN_EXP - e2; ndig -= DBL_MIN_EXP - e2;
# endif # endif
e2 -= ndig; e2 -= ndig;
mant = dbl2big(ldexp(dman, ndig)); mant = dbl2big(ldexp(dman, ndig));
point = ceil(e2*llog2); point = ceil(e2*llog2);
/* if (scm_verbose > 1) printf("mantissa = %g -> #x%s; e2 = %d -> %d; point = %d; ndig = %d -> %d\n", dman, CHARS(number2string(mant, MAKINUM(16))), e2+ndig, e2, point, dbl_mant_dig, ndig); */ /* if (scm_verbose > 1) printf("mantissa = %g -> #x%s; e2 = %d -> %d; point = %d; ndig = %d -> %d\n", dman, CHARS(number2string(mant, MAKINUM(16))), e2+ndig, e2, point, dbl_mant_dig, ndig); */
if (e2 >= 0) { if (e2 >= 0) {
/* try first with starved precision */ /* try first with starved precision */
{ {
num = scm_ash(mant, MAKINUM(e2 - point)); num = scm_ash(mant, MAKINUM(e2 - point));
bigrecy(mant);
quo = scm_round_quotient(num, VELTS(pows5)[(long) point]); quo = scm_round_quotient(num, VELTS(pows5)[(long) point]);
if (pmantexp2dbl(quo, point) != f) { if (pmantexp2dbl(quo, point) != f) {
bigrecy(quo); quo = num; int po2 = MAKINUM(1)==scm_logcount(mant);
if (quo != num) { bigrecy(quo); quo = num; }
num = scm_ash(quo, MAKINUM(1L)); num = scm_ash(quo, MAKINUM(1L));
bigrecy(quo); if (quo != num) bigrecy(quo);
quo = scm_round_quotient(num, VELTS(pows5)[(long) --point]); quo = scm_round_quotient(num, VELTS(pows5)[(long) --point]);
if (po2 && pmantexp2dbl(quo,point) != f) {
SCM quo1 = sum(quo, MAKINUM(1L));
if (quo != quo1) { bigrecy(quo); quo = quo1; }
if (pmantexp2dbl(quo,point) != f) {
if (quo != num) { bigrecy(quo); quo = num; };
num = scm_ash(quo,MAKINUM(1L));
quo = scm_round_quotient(num, VELTS(pows5)[(long) --point]);
}
}
} }
if (num != quo) bigrecy(num);
} }
} else { /* e2 <= 0 */ } else { /* e2 <= 0 */
/* try first with starved precision */ /* try first with starved precision */
{ {
SCM den = scm_ash(MAKINUM(1L), MAKINUM(point - e2)); SCM den = scm_ash(MAKINUM(1L), MAKINUM(point - e2));
num = product(mant, VELTS(pows5)[- (long) point]); num = product(mant, VELTS(pows5)[- (long) point]);
bigrecy(mant);
quo = scm_round_quotient(num, den); quo = scm_round_quotient(num, den);
if (pmantexp2dbl(quo, point) != f) { if (pmantexp2dbl(quo, point) != f) {
bigrecy(quo); quo = num; int po2 = MAKINUM(1L)==scm_logcount(mant);
if (quo != num) { bigrecy(quo); quo = num; }
point--; point--;
num = product(quo, MAKINUM(10)); num = product(quo, MAKINUM(10));
if (mant != MAKINUM(1)) bigrecy(quo); /* if (quo != num) bigrecy(quo); */
quo = scm_round_quotient(num, den); quo = scm_round_quotient(num, den);
} if (po2 && pmantexp2dbl(quo,point) != f) {
SCM quo1 = sum(quo, MAKINUM(1L));
if (quo != quo1) { bigrecy(quo); quo = quo1; }
if (pmantexp2dbl(quo,point) != f) {
if (quo != num) { bigrecy(quo); quo = num; }
point--;
num = product(quo, MAKINUM(10));
if (quo != num) bigrecy(quo);
quo = scm_round_quotient(num, den);
}
}
} else if (quo != num) bigrecy(num);
bigrecy(den); bigrecy(den);
} }
} }
bigrecy(num); if (mant != quo) bigrecy(mant);
a[ch++] = '.'; a[ch++] = '.';
/* if (sizeof(UBIGLONG)>=sizeof(double)) /\* Is ulong larger than mantissa? *\ / */ /* if (sizeof(UBIGLONG)>=sizeof(double)) /\* Is ulong larger than mantissa? *\ / */
/* ch += iulong2str(num2ulong(quo, (char *)ARG1, s_number2string), 10, &a[ch ]); */ /* ch += iulong2str(num2ulong(quo, (char *)ARG1, s_number2string), 10, &a[ch ]); */
/* else */ { /* else */ {
SCM str = number2string(quo, MAKINUM(10)); SCM str = number2string(quo, MAKINUM(10));
int len = LENGTH(str), i = 0; int len = LENGTH(str), i = 0;
bigrecy(quo); bigrecy(quo);
point += len - 1; point += len - 1;
while (i < len) a[ch++] = CHARS(str)[i++]; while (i < len) a[ch++] = CHARS(str)[i++];
strrecy(str); strrecy(str);
skipping to change at line 662 skipping to change at line 684
bigrecy(bmant); bigrecy(bmant);
return val; return val;
} }
double pmantexp2dbl(bmant, point) double pmantexp2dbl(bmant, point)
SCM bmant; SCM bmant;
int point; int point;
{ {
double ans; double ans;
if (BOOL_F == bmant) return REAL(scm_narn); if (BOOL_F == bmant) return REAL(scm_narn);
if (point >= 0) { if (point >= 0) {
SCM quo, num = product(bmant, VELTS(pows5)[(long) point]); if (point < 23 && INUM(scm_intlength(bmant)) <= dbl_mant_dig)
int bex = INUM(scm_intlength(num)) - dbl_mant_dig; return ldexp(num2dbl(bmant,ARG1,s_str2number) * dpows5[point], point);
if (bex > 0) { {
SCM den = scm_ash(MAKINUM(1L), MAKINUM(bex)); SCM quo, num = product(bmant, VELTS(pows5)[(long) point]);
quo = scm_round_quotient(num, den); int bex = INUM(scm_intlength(num)) - dbl_mant_dig;
bigrecy(den); if (bex > 0) {
} SCM den = scm_ash(MAKINUM(1L), MAKINUM(bex));
else quo = scm_round_quotient(num, den);
quo = scm_ash(num, MAKINUM(- bex)); bigrecy(den);
/* quo may not be a bignum */ }
if (INUMP(quo)) ans = ldexp((double)(INUM(quo)), bex + point); else
else { quo = scm_ash(num, MAKINUM(- bex));
sizet i = NUMDIGS(quo); /* quo may not be a bignum */
sizet j = i - (dbl_mant_dig + BITSPERDIG - 1)/BITSPERDIG; if (INUMP(quo)) ans = ldexp((double)(INUM(quo)), bex + point);
BIGDIG *digits = BDIGITS(quo); else {
ans = 0.0; sizet i = NUMDIGS(quo);
while (i-- > j) ans = digits[i] + ldexp(ans, BITSPERDIG); sizet j = (dbl_mant_dig + BITSPERDIG - 1)/BITSPERDIG;
bex += j * BITSPERDIG; BIGDIG *digits = BDIGITS(quo);
ans = ldexp(ans, bex + point); if (j < i) j = i - j;
else j = 0;
ans = 0.0;
while (i-- > j) ans = digits[i] + ldexp(ans, BITSPERDIG);
bex += j * BITSPERDIG;
ans = ldexp(ans, bex + point);
}
if (num != quo) bigrecy(quo);
if ((num != bmant) && (bmant != MAKINUM(1L))) bigrecy(num);
return ans;
} }
if (num != quo) bigrecy(quo); }
if (bmant != MAKINUM(1L)) bigrecy(num); if (-point < 23 && INUM(scm_intlength(bmant)) <= dbl_mant_dig)
return ans; return ldexp(num2dbl(bmant,ARG1,s_str2number) / dpows5[-point], point);
} else { {
int maxpow = LENGTH(pows5) - 1; int maxpow = LENGTH(pows5) - 1;
SCM scl = (-point <= maxpow) ? SCM num, quo, scl = (-point <= maxpow) ?
VELTS(pows5)[(long) -point] : VELTS(pows5)[(long) -point] :
product(VELTS(pows5)[(long)maxpow], VELTS(pows5)[(long)-point-maxpow]); product(VELTS(pows5)[(long)maxpow], VELTS(pows5)[(long)-point-maxpow]);
int mantlen = dbl_mant_dig;
int bex = /* bex < 0 */ int bex = /* bex < 0 */
INUM(scm_intlength(bmant)) - INUM(scm_intlength(scl)) - dbl_mant_dig; INUM(scm_intlength(bmant)) - INUM(scm_intlength(scl)) - mantlen;
SCM num = scm_ash(bmant, MAKINUM(-bex)); int tmp = bex + point + 1021 + mantlen;
SCM quo = scm_round_quotient(num, scl); if (tmp < 0) {
if (INUM(scm_intlength(quo)) > dbl_mant_dig) { bex -= tmp + 1;
mantlen += tmp;
}
num = scm_ash(bmant, MAKINUM(-bex));
quo = scm_round_quotient(num, scl);
if (INUM(scm_intlength(quo)) > mantlen) {
bex++; /* too many bits of quotient */ bex++; /* too many bits of quotient */
quo = scm_round_quotient(num, scm_ash(scl, MAKINUM(1L))); quo = scm_round_quotient(num, scm_ash(scl, MAKINUM(1L)));
} }
if (-point > maxpow) bigrecy(scl); if (-point > maxpow) bigrecy(scl);
bigrecy(num); if (num != quo) bigrecy(num);
ans = ldexp(int2dbl(quo), bex + point); ans = ldexp(int2dbl(quo), bex + point);
bigrecy(quo); bigrecy(quo);
return ans; return ans;
} }
} }
# else /* def BIGDIG */ # else /* def BIGDIG */
double mantexp2dbl(manstr, point) double mantexp2dbl(manstr, point)
SCM manstr; SCM manstr;
skipping to change at line 2479 skipping to change at line 2516
gencase: gencase:
{ {
SCM tz1 = z1; SCM tz1 = z1;
while(1) { while(1) {
SCM tmp; SCM tmp;
if (0L==iz2) break; if (0L==iz2) break;
if (1L==iz2) { if (1L==iz2) {
tmp = acc; tmp = acc;
acc = product(tmp, tz1); acc = product(tmp, tz1);
#ifdef BIGDIG #ifdef BIGDIG
bigrecy(tmp); if (acc != tmp) bigrecy(tmp);
#endif #endif
break; break;
} }
if (iz2 & 1) { if (iz2 & 1) {
tmp = acc; tmp = acc;
acc = product(tmp, tz1); acc = product(tmp, tz1);
#ifdef BIGDIG #ifdef BIGDIG
bigrecy(tmp); if (acc != tmp) bigrecy(tmp);
#endif #endif
} }
tmp = tz1; tmp = tz1;
tz1 = product(tmp, tmp); tz1 = product(tmp, tmp);
#ifdef BIGDIG #ifdef BIGDIG
if (tmp != z1 && tmp != acc) bigrecy(tmp); if (tmp != z1 && tmp != acc) bigrecy(tmp);
#endif #endif
iz2 >>= 1; iz2 >>= 1;
} }
#ifdef BIGDIG #ifdef BIGDIG
skipping to change at line 2834 skipping to change at line 2871
if (bex > 0) { if (bex > 0) {
SCM den = scm_ash(MAKINUM(1L), MAKINUM(bex)); SCM den = scm_ash(MAKINUM(1L), MAKINUM(bex));
quo = scm_round_quotient(num, den); quo = scm_round_quotient(num, den);
bigrecy(den); bigrecy(den);
} }
/* quo may not be a bignum */ /* quo may not be a bignum */
if (INUMP(quo)) if (INUMP(quo))
ans = ldexp((double)(INUM(quo)), bex); ans = ldexp((double)(INUM(quo)), bex);
else { else {
sizet i = NUMDIGS(quo); sizet i = NUMDIGS(quo);
sizet j = i - (dbl_mant_dig + BITSPERDIG - 1)/BITSPERDIG; sizet j = (dbl_mant_dig + BITSPERDIG - 1)/BITSPERDIG;
BIGDIG *digits = BDIGITS(quo); BIGDIG *digits = BDIGITS(quo);
if (j < 0) j = 0; if (j < i) j = i - j;
else j = 0;
while (i-- > j) ans = digits[i] + ldexp(ans, BITSPERDIG); while (i-- > j) ans = digits[i] + ldexp(ans, BITSPERDIG);
bex += j * BITSPERDIG; bex += j * BITSPERDIG;
if (bex > 0) ans = ldexp(ans, bex); if (bex > 0) ans = ldexp(ans, bex);
} }
if (quo != num) bigrecy(quo);
if (num != b) bigrecy(num); if (num != b) bigrecy(num);
if (quo != b) bigrecy(quo);
if (tc16_bigneg==TYP16(b)) return -ans; if (tc16_bigneg==TYP16(b)) return -ans;
return ans; return ans;
} }
} }
static SCM bigdblop(op, b, re, im) static SCM bigdblop(op, b, re, im)
int op; int op;
SCM b; SCM b;
double re, im; double re, im;
{ {
double bm = 0.0; double bm = 0.0;
skipping to change at line 3215 skipping to change at line 3253
double f = 0.1; double f = 0.1;
volatile double fsum = 1.0+f; volatile double fsum = 1.0+f;
while (fsum != 1.0) { while (fsum != 1.0) {
f /= 10.0; f /= 10.0;
if (++dblprec > 20) break; if (++dblprec > 20) break;
safe_add_1(f, &fsum); safe_add_1(f, &fsum);
} }
dblprec = dblprec-1; dblprec = dblprec-1;
} }
# endif /* DBL_DIG */ # endif /* DBL_DIG */
# else /* !BIGDIG */
{
int idx;
dpows5[0] = 1.0;
for (idx = 1; idx < 23; idx++) {
dpows5[idx] = 5*dpows5[idx-1];
}
}
# endif /* !BIGDIG */ # endif /* !BIGDIG */
# ifdef DBL_MANT_DIG # ifdef DBL_MANT_DIG
dbl_mant_dig = DBL_MANT_DIG; dbl_mant_dig = DBL_MANT_DIG;
# else # else
{ /* means we #defined it. */ { /* means we #defined it. */
volatile double fsum = 0.0; volatile double fsum = 0.0;
double eps = 1.0; double eps = 1.0;
int i = 0; int i = 0;
while (fsum != 1.0) { while (fsum != 1.0) {
eps = 0.5 * eps; eps = 0.5 * eps;
 End of changes. 25 change blocks. 
43 lines changed or deleted 89 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)