home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Usenet 1994 October
/
usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso
/
unix
/
volume27
/
calc-2.9.0
/
part11
< prev
next >
Wrap
Text File
|
1993-12-07
|
59KB
|
2,448 lines
Newsgroups: comp.sources.unix
From: dbell@canb.auug.org.au (David I. Bell)
Subject: v27i138: calc-2.9.0 - arbitrary precision C-like programmable calculator, Part11/19
References: <1.755316719.21314@gw.home.vix.com>
Sender: unix-sources-moderator@gw.home.vix.com
Approved: vixie@gw.home.vix.com
Submitted-By: dbell@canb.auug.org.au (David I. Bell)
Posting-Number: Volume 27, Issue 138
Archive-Name: calc-2.9.0/part11
#!/bin/sh
# this is part 11 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc2.9.0/value.c continued
#
CurArch=11
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
exit 1; fi
( read Scheck
if test "$Scheck" != $CurArch
then echo "Please unpack part $Scheck next!"
exit 1;
else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file calc2.9.0/value.c"
sed 's/^X//' << 'SHAR_EOF' >> calc2.9.0/value.c
X * Raise one value to another value's power, within the specified error.
X * Result is placed in the indicated location.
X */
Xvoid
Xpowervalue(v1, v2, v3, vres)
X VALUE *v1, *v2, *v3, *vres;
X{
X NUMBER *epsilon;
X COMPLEX *c, ctmp;
X
X vres->v_type = V_NULL;
X if (v3->v_type != V_NUM)
X math_error("Non-real epsilon value for power");
X epsilon = v3->v_num;
X if (qisneg(epsilon) || qiszero(epsilon))
X math_error("Non-positive epsilon value for power");
X switch (TWOVAL(v1->v_type, v2->v_type)) {
X case TWOVAL(V_NUM, V_NUM):
X vres->v_num = qpower(v1->v_num, v2->v_num, epsilon);
X vres->v_type = V_NUM;
X return;
X case TWOVAL(V_NUM, V_COM):
X ctmp.real = v1->v_num;
X ctmp.imag = &_qzero_;
X ctmp.links = 1;
X vres->v_com = cpower(&ctmp, v2->v_com, epsilon);
X break;
X case TWOVAL(V_COM, V_NUM):
X ctmp.real = v2->v_num;
X ctmp.imag = &_qzero_;
X ctmp.links = 1;
X vres->v_com = cpower(v1->v_com, &ctmp, epsilon);
X break;
X case TWOVAL(V_COM, V_COM):
X vres->v_com = cpower(v1->v_com, v2->v_com, epsilon);
X break;
X default:
X math_error("Illegal value for raising to power");
X }
X /*
X * Here for any complex result.
X */
X vres->v_type = V_COM;
X c = vres->v_com;
X if (!cisreal(c))
X return;
X vres->v_num = qlink(c->real);
X vres->v_type = V_NUM;
X comfree(c);
X}
X
X
X/*
X * Divide one arbitrary value by another one.
X * Result is placed in the indicated location.
X */
Xvoid
Xdivvalue(v1, v2, vres)
X VALUE *v1, *v2, *vres;
X{
X COMPLEX *c;
X COMPLEX ctmp;
X VALUE tmpval;
X
X vres->v_type = V_NULL;
X switch (TWOVAL(v1->v_type, v2->v_type)) {
X case TWOVAL(V_NUM, V_NUM):
X vres->v_num = qdiv(v1->v_num, v2->v_num);
X vres->v_type = V_NUM;
X return;
X case TWOVAL(V_COM, V_NUM):
X vres->v_com = cdivq(v1->v_com, v2->v_num);
X vres->v_type = V_COM;
X return;
X case TWOVAL(V_NUM, V_COM):
X if (qiszero(v1->v_num)) {
X vres->v_num = qlink(&_qzero_);
X vres->v_type = V_NUM;
X return;
X }
X ctmp.real = v1->v_num;
X ctmp.imag = &_qzero_;
X ctmp.links = 1;
X vres->v_com = cdiv(&ctmp, v2->v_com);
X vres->v_type = V_COM;
X return;
X case TWOVAL(V_COM, V_COM):
X vres->v_com = cdiv(v1->v_com, v2->v_com);
X vres->v_type = V_COM;
X c = vres->v_com;
X if (cisreal(c)) {
X vres->v_num = qlink(c->real);
X vres->v_type = V_NUM;
X comfree(c);
X }
X return;
X case TWOVAL(V_MAT, V_NUM):
X case TWOVAL(V_MAT, V_COM):
X invertvalue(v2, &tmpval);
X vres->v_mat = matmulval(v1->v_mat, &tmpval);
X vres->v_type = V_MAT;
X freevalue(&tmpval);
X return;
X default:
X if ((v1->v_type != V_OBJ) && (v2->v_type != V_OBJ))
X math_error("Non-compatible values for divide");
X *vres = objcall(OBJ_DIV, v1, v2, NULL_VALUE);
X return;
X }
X}
X
X
X/*
X * Divide one arbitrary value by another one keeping only the integer part.
X * Result is placed in the indicated location.
X */
Xvoid
Xquovalue(v1, v2, vres)
X VALUE *v1, *v2, *vres;
X{
X COMPLEX *c;
X
X vres->v_type = V_NULL;
X switch (TWOVAL(v1->v_type, v2->v_type)) {
X case TWOVAL(V_NUM, V_NUM):
X vres->v_num = qquo(v1->v_num, v2->v_num);
X vres->v_type = V_NUM;
X return;
X case TWOVAL(V_COM, V_NUM):
X vres->v_com = cquoq(v1->v_com, v2->v_num);
X vres->v_type = V_COM;
X c = vres->v_com;
X if (cisreal(c)) {
X vres->v_num = qlink(c->real);
X vres->v_type = V_NUM;
X comfree(c);
X }
X return;
X case TWOVAL(V_MAT, V_NUM):
X case TWOVAL(V_MAT, V_COM):
X vres->v_mat = matquoval(v1->v_mat, v2);
X vres->v_type = V_MAT;
X return;
X default:
X if ((v1->v_type != V_OBJ) && (v2->v_type != V_OBJ))
X math_error("Non-compatible values for quotient");
X *vres = objcall(OBJ_QUO, v1, v2, NULL_VALUE);
X return;
X }
X}
X
X
X/*
X * Divide one arbitrary value by another one keeping only the remainder.
X * Result is placed in the indicated location.
X */
Xvoid
Xmodvalue(v1, v2, vres)
X VALUE *v1, *v2, *vres;
X{
X COMPLEX *c;
X
X vres->v_type = V_NULL;
X switch (TWOVAL(v1->v_type, v2->v_type)) {
X case TWOVAL(V_NUM, V_NUM):
X vres->v_num = qmod(v1->v_num, v2->v_num);
X vres->v_type = V_NUM;
X return;
X case TWOVAL(V_COM, V_NUM):
X vres->v_com = cmodq(v1->v_com, v2->v_num);
X vres->v_type = V_COM;
X c = vres->v_com;
X if (cisreal(c)) {
X vres->v_num = qlink(c->real);
X vres->v_type = V_NUM;
X comfree(c);
X }
X return;
X case TWOVAL(V_MAT, V_NUM):
X case TWOVAL(V_MAT, V_COM):
X vres->v_mat = matmodval(v1->v_mat, v2);
X vres->v_type = V_MAT;
X return;
X default:
X if ((v1->v_type != V_OBJ) && (v2->v_type != V_OBJ))
X math_error("Non-compatible values for mod");
X *vres = objcall(OBJ_MOD, v1, v2, NULL_VALUE);
X return;
X }
X}
X
X
X/*
X * Test an arbitrary value to see if it is equal to "zero".
X * The definition of zero varies depending on the value type. For example,
X * the null string is "zero", and a matrix with zero values is "zero".
X * Returns TRUE if value is not equal to zero.
X */
XBOOL
Xtestvalue(vp)
X VALUE *vp;
X{
X VALUE val;
X
X switch (vp->v_type) {
X case V_NUM:
X return !qiszero(vp->v_num);
X case V_COM:
X return !ciszero(vp->v_com);
X case V_STR:
X return (vp->v_str[0] != '\0');
X case V_MAT:
X return mattest(vp->v_mat);
X case V_LIST:
X return (vp->v_list->l_count != 0);
X case V_ASSOC:
X return (vp->v_assoc->a_count != 0);
X case V_FILE:
X return validid(vp->v_file);
X case V_NULL:
X return FALSE;
X case V_OBJ:
X val = objcall(OBJ_TEST, vp, NULL_VALUE, NULL_VALUE);
X return (val.v_int != 0);
X default:
X return TRUE;
X }
X}
X
X
X/*
X * Compare two values for equality.
X * Returns TRUE if the two values differ.
X */
XBOOL
Xcomparevalue(v1, v2)
X VALUE *v1, *v2;
X{
X int r;
X VALUE val;
X
X if ((v1->v_type == V_OBJ) || (v2->v_type == V_OBJ)) {
X val = objcall(OBJ_CMP, v1, v2, NULL_VALUE);
X return (val.v_int != 0);
X }
X if (v1 == v2)
X return FALSE;
X if (v1->v_type != v2->v_type)
X return TRUE;
X switch (v1->v_type) {
X case V_NUM:
X r = qcmp(v1->v_num, v2->v_num);
X break;
X case V_COM:
X r = ccmp(v1->v_com, v2->v_com);
X break;
X case V_STR:
X r = ((v1->v_str != v2->v_str) &&
X ((v1->v_str[0] - v2->v_str[0]) ||
X strcmp(v1->v_str, v2->v_str)));
X break;
X case V_MAT:
X r = matcmp(v1->v_mat, v2->v_mat);
X break;
X case V_LIST:
X r = listcmp(v1->v_list, v2->v_list);
X break;
X case V_ASSOC:
X r = assoccmp(v1->v_assoc, v2->v_assoc);
X break;
X case V_NULL:
X r = FALSE;
X break;
X case V_FILE:
X r = (v1->v_file != v2->v_file);
X break;
X default:
X math_error("Illegal values for comparevalue");
X }
X return (r != 0);
X}
X
X
X/*
X * Compare two values for their relative values.
X * Returns minus one if the first value is less than the second one,
X * one if the first value is greater than the second one, and
X * zero if they are equal.
X */
XFLAG
Xrelvalue(v1, v2)
X VALUE *v1, *v2;
X{
X int r;
X VALUE val;
X
X if ((v1->v_type == V_OBJ) || (v2->v_type == V_OBJ)) {
X val = objcall(OBJ_REL, v1, v2, NULL_VALUE);
X return val.v_int;
X }
X if (v1 == v2)
X return 0;
X if (v1->v_type != v2->v_type)
X math_error("Relative comparison of differing types");
X switch (v1->v_type) {
X case V_NUM:
X r = qrel(v1->v_num, v2->v_num);
X break;
X case V_STR:
X r = strcmp(v1->v_str, v2->v_str);
X break;
X case V_NULL:
X r = 0;
X break;
X default:
X math_error("Illegal value for relative comparison");
X }
X if (r < 0)
X return -1;
X return (r != 0);
X}
X
X
X/*
X * Calculate a hash value for a value.
X * The hash does not have to be a perfect one, it is only used for
X * making associations faster.
X */
XHASH
Xhashvalue(vp)
X VALUE *vp;
X{
X switch (vp->v_type) {
X case V_INT:
X return ((long) vp->v_int);
X break;
X case V_NUM:
X return qhash(vp->v_num);
X break;
X case V_COM:
X return chash(vp->v_com);
X break;
X case V_STR:
X return hashstr(vp->v_str);
X break;
X case V_NULL:
X return 0;
X break;
X case V_OBJ:
X return objhash(vp->v_obj);
X break;
X case V_LIST:
X return listhash(vp->v_list);
X break;
X case V_ASSOC:
X return assochash(vp->v_assoc);
X break;
X case V_MAT:
X return mathash(vp->v_mat);
X break;
X case V_FILE:
X return ((long) vp->v_file);
X break;
X default:
X math_error("Hashing unknown value");
X }
X return 0;
X}
X
X
X/*
X * Print the value of a descriptor in one of several formats.
X * If flags contains PRINT_SHORT, then elements of arrays and lists
X * will not be printed. If flags contains PRINT_UNAMBIG, then quotes
X * are placed around strings and the null value is explicitly printed.
X */
Xvoid
Xprintvalue(vp, flags)
X VALUE *vp;
X{
X switch (vp->v_type) {
X case V_NUM:
X qprintnum(vp->v_num, MODE_DEFAULT);
X break;
X case V_COM:
X comprint(vp->v_com);
X break;
X case V_STR:
X if (flags & PRINT_UNAMBIG)
X math_chr('\"');
X math_str(vp->v_str);
X if (flags & PRINT_UNAMBIG)
X math_chr('\"');
X break;
X case V_NULL:
X if (flags & PRINT_UNAMBIG)
X math_str("NULL");
X break;
X case V_OBJ:
X (void) objcall(OBJ_PRINT, vp, NULL_VALUE, NULL_VALUE);
X break;
X case V_LIST:
X listprint(vp->v_list,
X ((flags & PRINT_SHORT) ? 0L : maxprint));
X break;
X case V_ASSOC:
X assocprint(vp->v_assoc,
X ((flags & PRINT_SHORT) ? 0L : maxprint));
X break;
X case V_MAT:
X matprint(vp->v_mat,
X ((flags & PRINT_SHORT) ? 0L : maxprint));
X break;
X case V_FILE:
X printid(vp->v_file, flags);
X break;
X default:
X math_error("Printing unknown value");
X }
X}
X
X/* END CODE */
SHAR_EOF
echo "File calc2.9.0/value.c is complete"
chmod 0644 calc2.9.0/value.c || echo "restore of calc2.9.0/value.c fails"
set `wc -c calc2.9.0/value.c`;Sum=$1
if test "$Sum" != "29836"
then echo original size 29836, current size $Sum;fi
echo "x - extracting calc2.9.0/value.h (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/value.h &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Definitions of general values and related routines used by the calculator.
X */
X
X#ifndef VALUE_H
X#define VALUE_H
X
X#include "cmath.h"
X
X
X#define MAXDIM 4 /* maximum number of dimensions in matrices */
X#define USUAL_ELEMENTS 4 /* usual number of elements for objects */
X
X
X/*
X * Flags to modify results from the printvalue routine.
X * These flags are OR'd together.
X */
X#define PRINT_NORMAL 0x00 /* print in normal manner */
X#define PRINT_SHORT 0x01 /* print in short format (no elements) */
X#define PRINT_UNAMBIG 0x02 /* print in non-ambiguous manner */
X
X
X/*
X * Definition of values of various types.
X */
Xtypedef struct value VALUE;
Xtypedef struct object OBJECT;
Xtypedef struct matrix MATRIX;
Xtypedef struct list LIST;
Xtypedef struct assoc ASSOC;
Xtypedef long FILEID;
X
X
Xstruct value {
X short v_type; /* type of value */
X short v_subtype; /* other data related to some types */
X union {
X long vv_int; /* small integer value */
X FILEID vv_file; /* id of opened file */
X NUMBER *vv_num; /* arbitrary sized numeric value */
X COMPLEX *vv_com; /* complex number */
X VALUE *vv_addr; /* address of variable value */
X MATRIX *vv_mat; /* address of matrix */
X LIST *vv_list; /* address of list */
X ASSOC *vv_assoc; /* address of association */
X OBJECT *vv_obj; /* address of object */
X char *vv_str; /* string value */
X } v_union;
X};
X
X
X/*
X * For ease in referencing
X */
X#define v_int v_union.vv_int
X#define v_file v_union.vv_file
X#define v_num v_union.vv_num
X#define v_com v_union.vv_com
X#define v_addr v_union.vv_addr
X#define v_str v_union.vv_str
X#define v_mat v_union.vv_mat
X#define v_list v_union.vv_list
X#define v_assoc v_union.vv_assoc
X#define v_obj v_union.vv_obj
X#define v_valid v_union.vv_int
X
X
X/*
X * Value types.
X */
X#define V_NULL 0 /* null value */
X#define V_INT 1 /* normal integer */
X#define V_NUM 2 /* number */
X#define V_COM 3 /* complex number */
X#define V_ADDR 4 /* address of variable value */
X#define V_STR 5 /* address of string */
X#define V_MAT 6 /* address of matrix structure */
X#define V_LIST 7 /* address of list structure */
X#define V_ASSOC 8 /* address of association structure */
X#define V_OBJ 9 /* address of object structure */
X#define V_FILE 10 /* opened file id */
X#define V_MAX 10 /* highest legal value */
X
X#define V_STRLITERAL 0 /* string subtype for literal str */
X#define V_STRALLOC 1 /* string subtype for allocated str */
X
X#define TWOVAL(a,b) ((a) * (V_MAX+1) + (b)) /* for switch of two values */
X
X#define NULL_VALUE ((VALUE *) 0)
X
X
Xextern void freevalue MATH_PROTO((VALUE *vp));
Xextern void copyvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void negvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void addvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void subvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void mulvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void squarevalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void invertvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void roundvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void broundvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void intvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void fracvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void incvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void decvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void conjvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void sqrtvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void rootvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *v3,
X VALUE *vres));
Xextern void absvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void normvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void shiftvalue MATH_PROTO((VALUE *v1, VALUE *v2, BOOL rightshift,
X VALUE *vres));
Xextern void scalevalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void powivalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void powervalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *v3,
X VALUE *vres));
Xextern void divvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void quovalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void modvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern BOOL testvalue MATH_PROTO((VALUE *vp));
Xextern BOOL comparevalue MATH_PROTO((VALUE *v1, VALUE *v2));
Xextern FLAG relvalue MATH_PROTO((VALUE *v1, VALUE *v2));
Xextern HASH hashvalue MATH_PROTO((VALUE *vp));
Xextern void printvalue MATH_PROTO((VALUE *vp, int flags));
X
X
X
X/*
X * Structure of a matrix.
X */
Xstruct matrix {
X long m_dim; /* dimension of matrix */
X long m_size; /* total number of elements */
X long m_min[MAXDIM]; /* minimum bound for indices */
X long m_max[MAXDIM]; /* maximum bound for indices */
X VALUE m_table[1]; /* actually varying length table */
X};
X
X#define matsize(n) (sizeof(MATRIX) - sizeof(VALUE) + ((n) * sizeof(VALUE)))
X
X
Xextern MATRIX *matadd MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern MATRIX *matsub MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern MATRIX *matmul MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern MATRIX *matneg MATH_PROTO((MATRIX *m));
Xextern MATRIX *matalloc MATH_PROTO((long size));
Xextern MATRIX *matcopy MATH_PROTO((MATRIX *m));
Xextern MATRIX *matsquare MATH_PROTO((MATRIX *m));
Xextern MATRIX *matinv MATH_PROTO((MATRIX *m));
Xextern MATRIX *matscale MATH_PROTO((MATRIX *m, long n));
Xextern MATRIX *matshift MATH_PROTO((MATRIX *m, long n));
Xextern MATRIX *matmulval MATH_PROTO((MATRIX *m, VALUE *vp));
Xextern MATRIX *matpowi MATH_PROTO((MATRIX *m, NUMBER *q));
Xextern MATRIX *matconj MATH_PROTO((MATRIX *m));
Xextern MATRIX *matquoval MATH_PROTO((MATRIX *m, VALUE *vp));
Xextern MATRIX *matmodval MATH_PROTO((MATRIX *m, VALUE *vp));
Xextern MATRIX *matint MATH_PROTO((MATRIX *m));
Xextern MATRIX *matfrac MATH_PROTO((MATRIX *m));
Xextern MATRIX *matround MATH_PROTO((MATRIX *m, long places));
Xextern MATRIX *matbround MATH_PROTO((MATRIX *m, long places));
Xextern MATRIX *mattrans MATH_PROTO((MATRIX *m));
Xextern MATRIX *matcross MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern BOOL mattest MATH_PROTO((MATRIX *m));
Xextern BOOL matcmp MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern long matsearch MATH_PROTO((MATRIX *m, VALUE *vp, long index));
Xextern long matrsearch MATH_PROTO((MATRIX *m, VALUE *vp, long index));
Xextern HASH mathash MATH_PROTO((MATRIX *m));
Xextern VALUE matdet MATH_PROTO((MATRIX *m));
Xextern VALUE matdot MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern void matfill MATH_PROTO((MATRIX *m, VALUE *v1, VALUE *v2));
Xextern void matfree MATH_PROTO((MATRIX *m));
Xextern void matprint MATH_PROTO((MATRIX *m, long max_print));
Xextern VALUE *matindex MATH_PROTO((MATRIX *mp, BOOL create, long dim,
X VALUE *indices));
X
X
X#if 0
Xextern BOOL matisident MATH_PROTO((MATRIX *m));
X#endif
X
X
X
X/*
X * List definitions.
X * An individual list element.
X */
Xtypedef struct listelem LISTELEM;
Xstruct listelem {
X LISTELEM *e_next; /* next element in list (or NULL) */
X LISTELEM *e_prev; /* previous element in list (or NULL) */
X VALUE e_value; /* value of this element */
X};
X
X
X/*
X * Structure for a list of elements.
X */
Xstruct list {
X LISTELEM *l_first; /* first list element (or NULL) */
X LISTELEM *l_last; /* last list element (or NULL) */
X LISTELEM *l_cache; /* cached list element (or NULL) */
X long l_cacheindex; /* index of cached element (or undefined) */
X long l_count; /* total number of elements in the list */
X};
X
X
Xextern void insertlistfirst MATH_PROTO((LIST *lp, VALUE *vp));
Xextern void insertlistlast MATH_PROTO((LIST *lp, VALUE *vp));
Xextern void insertlistmiddle MATH_PROTO((LIST *lp, long index, VALUE *vp));
Xextern void removelistfirst MATH_PROTO((LIST *lp, VALUE *vp));
Xextern void removelistlast MATH_PROTO((LIST *lp, VALUE *vp));
Xextern void removelistmiddle MATH_PROTO((LIST *lp, long index, VALUE *vp));
Xextern void listfree MATH_PROTO((LIST *lp));
Xextern void listprint MATH_PROTO((LIST *lp, long max_print));
Xextern long listsearch MATH_PROTO((LIST *lp, VALUE *vp, long index));
Xextern long listrsearch MATH_PROTO((LIST *lp, VALUE *vp, long index));
Xextern HASH listhash MATH_PROTO((LIST *lp));
Xextern BOOL listcmp MATH_PROTO((LIST *lp1, LIST *lp2));
Xextern VALUE *listfindex MATH_PROTO((LIST *lp, long index));
Xextern LIST *listalloc MATH_PROTO((void));
Xextern LIST *listcopy MATH_PROTO((LIST *lp));
X
X
X
X/*
X * Structures for associations.
X * Associations are "indexed" by one or more arbitrary values, and are
X * stored in a hash table with their hash values for quick indexing.
X */
Xtypedef struct assocelem ASSOCELEM;
Xstruct assocelem {
X ASSOCELEM *e_next; /* next element in list (or NULL) */
X long e_dim; /* dimension of indexing for this element */
X HASH e_hash; /* hash value for this element */
X VALUE e_value; /* value of association */
X VALUE e_indices[1]; /* index values (variable length) */
X};
X
X
Xstruct assoc {
X long a_count; /* number of elements in the association */
X long a_size; /* current size of association hash table */
X ASSOCELEM **a_table; /* current hash table for elements */
X};
X
X
Xextern ASSOC *assocalloc MATH_PROTO((long initsize));
Xextern ASSOC *assoccopy MATH_PROTO((ASSOC *ap));
Xextern void assocfree MATH_PROTO((ASSOC *ap));
Xextern void assocprint MATH_PROTO((ASSOC *ap, long max_print));
Xextern long assocsearch MATH_PROTO((ASSOC *ap, VALUE *vp, long index));
Xextern long assocrsearch MATH_PROTO((ASSOC *ap, VALUE *vp, long index));
Xextern HASH assochash MATH_PROTO((ASSOC *ap));
Xextern BOOL assoccmp MATH_PROTO((ASSOC *ap1, ASSOC *ap2));
Xextern VALUE *assocfindex MATH_PROTO((ASSOC *ap, long index));
Xextern VALUE *associndex MATH_PROTO((ASSOC *ap, BOOL create, long dim,
X VALUE *indices));
X
X
X/*
X * Object actions.
X */
X#define OBJ_PRINT 0 /* print the value */
X#define OBJ_ONE 1 /* create the multiplicative identity */
X#define OBJ_TEST 2 /* test a value for "zero" */
X#define OBJ_ADD 3 /* add two values */
X#define OBJ_SUB 4 /* subtrace one value from another */
X#define OBJ_NEG 5 /* negate a value */
X#define OBJ_MUL 6 /* multiply two values */
X#define OBJ_DIV 7 /* divide one value by another */
X#define OBJ_INV 8 /* invert a value */
X#define OBJ_ABS 9 /* take absolute value of value */
X#define OBJ_NORM 10 /* take the norm of a value */
X#define OBJ_CONJ 11 /* take the conjugate of a value */
X#define OBJ_POW 12 /* take the power function */
X#define OBJ_SGN 13 /* return the sign of a value */
X#define OBJ_CMP 14 /* compare two values for equality */
X#define OBJ_REL 15 /* compare two values for inequality */
X#define OBJ_QUO 16 /* integer quotient of values */
X#define OBJ_MOD 17 /* remainder of division of values */
X#define OBJ_INT 18 /* integer part of */
X#define OBJ_FRAC 19 /* fractional part of */
X#define OBJ_INC 20 /* increment by one */
X#define OBJ_DEC 21 /* decrement by one */
X#define OBJ_SQUARE 22 /* square value */
X#define OBJ_SCALE 23 /* scale by power of two */
X#define OBJ_SHIFT 24 /* shift left (or right) by number of bits */
X#define OBJ_ROUND 25 /* round to specified decimal places */
X#define OBJ_BROUND 26 /* round to specified binary places */
X#define OBJ_ROOT 27 /* take nth root of value */
X#define OBJ_SQRT 28 /* take square root of value */
X#define OBJ_MAXFUNC 28 /* highest function */
X
X
X/*
X * Definition of an object type.
X * This is actually a varying sized structure.
X */
Xtypedef struct {
X char *name; /* name of object */
X int count; /* number of elements defined */
X int actions[OBJ_MAXFUNC+1]; /* function indices for actions */
X int elements[1]; /* element indexes (MUST BE LAST) */
X} OBJECTACTIONS;
X
X#define objectactionsize(elements) \
X (sizeof(OBJECTACTIONS) + ((elements) - 1) * sizeof(int))
X
X
X/*
X * Structure of an object.
X * This is actually a varying sized structure.
X * However, there are always at least USUAL_ELEMENTS values in the object.
X */
Xstruct object {
X OBJECTACTIONS *o_actions; /* action table for this object */
X VALUE o_table[USUAL_ELEMENTS]; /* object values (MUST BE LAST) */
X};
X
X#define objectsize(elements) \
X (sizeof(OBJECT) + ((elements) - USUAL_ELEMENTS) * sizeof(VALUE))
X
X
Xextern OBJECT *objcopy MATH_PROTO((OBJECT *op));
Xextern OBJECT *objalloc MATH_PROTO((long index));
Xextern VALUE objcall MATH_PROTO((int action, VALUE *v1, VALUE *v2, VALUE *v3));
Xextern void objfree MATH_PROTO((OBJECT *op));
Xextern void objuncache MATH_PROTO((void));
Xextern int addelement MATH_PROTO((char *name));
Xextern void defineobject MATH_PROTO((char *name, int indices[], int count));
Xextern int checkobject MATH_PROTO((char *name));
Xextern void showobjfuncs MATH_PROTO((void));
Xextern int findelement MATH_PROTO((char *name));
Xextern int objoffset MATH_PROTO((OBJECT *op, long index));
Xextern HASH objhash MATH_PROTO((OBJECT *op));
X
X#endif
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/value.h || echo "restore of calc2.9.0/value.h fails"
set `wc -c calc2.9.0/value.h`;Sum=$1
if test "$Sum" != "12825"
then echo original size 12825, current size $Sum;fi
echo "x - extracting calc2.9.0/version.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/version.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * version - determine the version of calc
X */
X
X#include "calc.h"
X
X#define MAJOR_VER 2 /* major version */
X#define MINOR_VER 9 /* minor version */
X#define PATCH_LEVEL 0 /* patch level */
X
X
Xvoid
Xversion(stream)
X FILE *stream; /* stream to write version on */
X{
X fprintf(stream,
X "C-style arbitrary precision calculator (version %d.%d.%d)\n",
X MAJOR_VER, MINOR_VER, PATCH_LEVEL);
X}
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/version.c || echo "restore of calc2.9.0/version.c fails"
set `wc -c calc2.9.0/version.c`;Sum=$1
if test "$Sum" != "564"
then echo original size 564, current size $Sum;fi
echo "x - extracting calc2.9.0/zfunc.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/zfunc.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Extended precision integral arithmetic non-primitive routines
X */
X
X#include "zmath.h"
X
Xstatic ZVALUE primeprod; /* product of primes under 100 */
XZVALUE _tenpowers_[32]; /* table of 10^2^n */
X
X
X/*
X * Compute the factorial of a number.
X */
Xvoid
Xzfact(z, dest)
X ZVALUE z, *dest;
X{
X long ptwo; /* count of powers of two */
X long n; /* current multiplication value */
X long m; /* reduced multiplication value */
X long mul; /* collected value to multiply by */
X ZVALUE res, temp;
X
X if (zisneg(z))
X math_error("Negative argument for factorial");
X if (zisbig(z))
X math_error("Very large factorial");
X n = (zistiny(z) ? z1tol(z) : z2tol(z));
X ptwo = 0;
X mul = 1;
X res = _one_;
X /*
X * Multiply numbers together, but squeeze out all powers of two.
X * We will put them back in at the end. Also collect multiple
X * numbers together until there is a risk of overflow.
X */
X for (; n > 1; n--) {
X for (m = n; ((m & 0x1) == 0); m >>= 1)
X ptwo++;
X mul *= m;
X if (mul < BASE1/2)
X continue;
X zmuli(res, mul, &temp);
X zfree(res);
X res = temp;
X mul = 1;
X }
X /*
X * Multiply by the remaining value, then scale result by
X * the proper power of two.
X */
X if (mul > 1) {
X zmuli(res, mul, &temp);
X zfree(res);
X res = temp;
X }
X zshift(res, ptwo, &temp);
X zfree(res);
X *dest = temp;
X}
X
X
X/*
X * Compute the product of the primes up to the specified number.
X */
Xvoid
Xzpfact(z, dest)
X ZVALUE z, *dest;
X{
X long n; /* limiting number to multiply by */
X long p; /* current prime */
X long i; /* test value */
X long mul; /* collected value to multiply by */
X ZVALUE res, temp;
X
X if (zisneg(z))
X math_error("Negative argument for factorial");
X if (zisbig(z))
X math_error("Very large factorial");
X n = (zistiny(z) ? z1tol(z) : z2tol(z));
X /*
X * Multiply by the primes in order, collecting multiple numbers
X * together until there is a chance of overflow.
X */
X mul = 1 + (n > 1);
X res = _one_;
X for (p = 3; p <= n; p += 2) {
X for (i = 3; (i * i) <= p; i += 2) {
X if ((p % i) == 0)
X goto next;
X }
X mul *= p;
X if (mul < BASE1/2)
X continue;
X zmuli(res, mul, &temp);
X zfree(res);
X res = temp;
X mul = 1;
Xnext: ;
X }
X /*
X * Multiply by the final value if any.
X */
X if (mul > 1) {
X zmuli(res, mul, &temp);
X zfree(res);
X res = temp;
X }
X *dest = res;
X}
X
X
X/*
X * Compute the least common multiple of all the numbers up to the
X * specified number.
X */
Xvoid
Xzlcmfact(z, dest)
X ZVALUE z, *dest;
X{
X long n; /* limiting number to multiply by */
X long p; /* current prime */
X long pp; /* power of prime */
X long i; /* test value */
X ZVALUE res, temp;
X
X if (zisneg(z) || ziszero(z))
X math_error("Non-positive argument for lcmfact");
X if (zisbig(z))
X math_error("Very large lcmfact");
X n = (zistiny(z) ? z1tol(z) : z2tol(z));
X /*
X * Multiply by powers of the necessary odd primes in order.
X * The power for each prime is the highest one which is not
X * more than the specified number.
X */
X res = _one_;
X for (p = 3; p <= n; p += 2) {
X for (i = 3; (i * i) <= p; i += 2) {
X if ((p % i) == 0)
X goto next;
X }
X i = p;
X while (i <= n) {
X pp = i;
X i *= p;
X }
X zmuli(res, pp, &temp);
X zfree(res);
X res = temp;
Xnext: ;
X }
X /*
X * Finish by scaling by the necessary power of two.
X */
X zshift(res, zhighbit(z), dest);
X zfree(res);
X}
X
X
X/*
X * Compute the permutation function M! / (M - N)!.
X */
Xvoid
Xzperm(z1, z2, res)
X ZVALUE z1, z2, *res;
X{
X long count;
X ZVALUE cur, tmp, ans;
X
X if (zisneg(z1) || zisneg(z2))
X math_error("Negative argument for permutation");
X if (zrel(z1, z2) < 0)
X math_error("Second arg larger than first in permutation");
X if (zisbig(z2))
X math_error("Very large permutation");
X count = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
X zcopy(z1, &ans);
X zsub(z1, _one_, &cur);
X while (--count > 0) {
X zmul(ans, cur, &tmp);
X zfree(ans);
X ans = tmp;
X zsub(cur, _one_, &tmp);
X zfree(cur);
X cur = tmp;
X }
X zfree(cur);
X *res = ans;
X}
X
X
X/*
X * Compute the combinatorial function M! / ( N! * (M - N)! ).
X */
Xvoid
Xzcomb(z1, z2, res)
X ZVALUE z1, z2, *res;
X{
X ZVALUE ans;
X ZVALUE mul, div, temp;
X FULL count, i;
X HALF dh[2];
X
X if (zisneg(z1) || zisneg(z2))
X math_error("Negative argument for combinatorial");
X zsub(z1, z2, &temp);
X if (zisneg(temp)) {
X zfree(temp);
X math_error("Second arg larger than first for combinatorial");
X }
X if (zisbig(z2) && zisbig(temp)) {
X zfree(temp);
X math_error("Very large combinatorial");
X }
X count = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
X i = (zistiny(temp) ? z1tol(temp) : z2tol(temp));
X if (zisbig(z2) || (!zisbig(temp) && (i < count)))
X count = i;
X zfree(temp);
X mul = z1;
X div.sign = 0;
X div.v = dh;
X ans = _one_;
X for (i = 1; i <= count; i++) {
X dh[0] = i & BASE1;
X dh[1] = i / BASE;
X div.len = 1 + (dh[1] != 0);
X zmul(ans, mul, &temp);
X zfree(ans);
X zquo(temp, div, &ans);
X zfree(temp);
X zsub(mul, _one_, &temp);
X if (mul.v != z1.v)
X zfree(mul);
X mul = temp;
X }
X if (mul.v != z1.v)
X zfree(mul);
X *res = ans;
X}
X
X
X/*
X * Perform a probabilistic primality test (algorithm P in Knuth).
X * Returns FALSE if definitely not prime, or TRUE if probably prime.
X * Count determines how many times to check for primality.
X * The chance of a non-prime passing this test is less than (1/4)^count.
X * For example, a count of 100 fails for only 1 in 10^60 numbers.
X */
XBOOL
Xzprimetest(z, count)
X ZVALUE z; /* number to test for primeness */
X long count;
X{
X long ij, ik, ix;
X ZVALUE zm1, z1, z2, z3, ztmp;
X HALF val[2];
X
X z.sign = 0;
X if (ziseven(z)) /* if even, not prime if not 2 */
X return (zistwo(z) != 0);
X /*
X * See if the number is small, and is either a small prime,
X * or is divisible by a small prime.
X */
X if (zistiny(z) && (*z.v <= (HALF)(101*101-1))) {
X ix = *z.v;
X for (ik = 3; (ik <= 97) && ((ik * ik) <= ix); ik += 2)
X if ((ix % ik) == 0)
X return FALSE;
X return TRUE;
X }
X /*
X * See if the number is divisible by one of the primes 3, 5,
X * 7, 11, or 13. This is a very easy check.
X */
X ij = zmodi(z, 15015L);
X if (!(ij % 3) || !(ij % 5) || !(ij % 7) || !(ij % 11) || !(ij % 13))
X return FALSE;
X /*
X * Check the gcd of the number and the product of more of the first
X * few odd primes. We must build the prime product on the first call.
X */
X ztmp.sign = 0;
X ztmp.len = 1;
X ztmp.v = val;
X if (primeprod.len == 0) {
X val[0] = 101;
X zpfact(ztmp, &primeprod);
X }
X zgcd(z, primeprod, &z1);
X if (!zisunit(z1)) {
X zfree(z1);
X return FALSE;
X }
X zfree(z1);
X /*
X * Not divisible by a small prime, so onward with the real test.
X * Make sure the count is limited by the number of odd numbers between
X * three and the number being tested.
X */
X ix = ((zistiny(z) ? z1tol(z) : z2tol(z) - 3) / 2);
X if (count > ix) count = ix;
X zsub(z, _one_, &zm1);
X ik = zlowbit(zm1);
X zshift(zm1, -ik, &z1);
X /*
X * Loop over various "random" numbers, testing each one.
X * These numbers are the odd numbers starting from three.
X */
X for (ix = 0; ix < count; ix++) {
X val[0] = (ix * 2) + 3;
X ij = 0;
X zpowermod(ztmp, z1, z, &z3);
X for (;;) {
X if (zisone(z3)) {
X if (ij) /* number is definitely not prime */
X goto notprime;
X break;
X }
X if (zcmp(z3, zm1) == 0)
X break;
X if (++ij >= ik)
X goto notprime; /* number is definitely not prime */
X zsquare(z3, &z2);
X zfree(z3);
X zmod(z2, z, &z3);
X zfree(z2);
X }
X zfree(z3);
X }
X zfree(zm1);
X zfree(z1);
X return TRUE; /* number might be prime */
X
Xnotprime:
X zfree(z3);
X zfree(zm1);
X zfree(z1);
X return FALSE;
X}
X
X
X/*
X * Compute the Jacobi function (p / q) for odd q.
X * If q is prime then the result is:
X * 1 if p == x^2 (mod q) for some x.
X * -1 otherwise.
X * If q is not prime, then the result is not meaningful if it is 1.
X * This function returns 0 if q is even or q < 0.
X */
XFLAG
Xzjacobi(z1, z2)
X ZVALUE z1, z2;
X{
X ZVALUE p, q, tmp;
X long lowbit;
X int val;
X
X if (ziseven(z2) || zisneg(z2))
X return 0;
X val = 1;
X if (ziszero(z1) || zisone(z1))
X return val;
X if (zisunit(z1)) {
X if ((*z2.v - 1) & 0x2)
X val = -val;
X return val;
X }
X zcopy(z1, &p);
X zcopy(z2, &q);
X for (;;) {
X zmod(p, q, &tmp);
X zfree(p);
X p = tmp;
X if (ziszero(p)) {
X zfree(p);
X p = _one_;
X }
X if (ziseven(p)) {
X lowbit = zlowbit(p);
X zshift(p, -lowbit, &tmp);
X zfree(p);
X p = tmp;
X if ((lowbit & 1) && (((*q.v & 0x7) == 3) || ((*q.v & 0x7) == 5)))
X val = -val;
X }
X if (zisunit(p)) {
X zfree(p);
X zfree(q);
X return val;
X }
X if ((*p.v & *q.v & 0x3) == 3)
X val = -val;
X tmp = q;
X q = p;
X p = tmp;
X }
X}
X
X
X/*
X * Return the Fibonacci number F(n).
X * This is evaluated by recursively using the formulas:
X * F(2N+1) = F(N+1)^2 + F(N)^2
X * and
X * F(2N) = F(N+1)^2 - F(N-1)^2
X */
Xvoid
Xzfib(z, res)
X ZVALUE z, *res;
X{
X unsigned long i;
X long n;
X int sign;
X ZVALUE fnm1, fn, fnp1; /* consecutive fibonacci values */
X ZVALUE t1, t2, t3;
X
X if (zisbig(z))
X math_error("Very large Fibonacci number");
X n = (zistiny(z) ? z1tol(z) : z2tol(z));
X if (n == 0) {
X *res = _zero_;
X return;
X }
X sign = z.sign && ((n & 0x1) == 0);
X if (n <= 2) {
X *res = _one_;
X res->sign = (BOOL)sign;
X return;
X }
X i = TOPFULL;
X while ((i & n) == 0)
X i >>= 1L;
X i >>= 1L;
X fnm1 = _zero_;
X fn = _one_;
X fnp1 = _one_;
X while (i) {
X zsquare(fnm1, &t1);
X zsquare(fn, &t2);
X zsquare(fnp1, &t3);
X zfree(fnm1);
X zfree(fn);
X zfree(fnp1);
X zadd(t2, t3, &fnp1);
X zsub(t3, t1, &fn);
X zfree(t1);
X zfree(t2);
X zfree(t3);
X if (i & n) {
X fnm1 = fn;
X fn = fnp1;
X zadd(fnm1, fn, &fnp1);
X } else
X zsub(fnp1, fn, &fnm1);
X i >>= 1L;
X }
X zfree(fnm1);
X zfree(fnp1);
X *res = fn;
X res->sign = (BOOL)sign;
X}
X
X
X/*
X * Compute the result of raising one number to the power of another
X * The second number is assumed to be non-negative.
X * It cannot be too large except for trivial cases.
X */
Xvoid
Xzpowi(z1, z2, res)
X ZVALUE z1, z2, *res;
X{
X int sign; /* final sign of number */
X unsigned long power; /* power to raise to */
X unsigned long bit; /* current bit value */
X long twos; /* count of times 2 is in result */
X ZVALUE ans, temp;
X
X sign = (z1.sign && zisodd(z2));
X z1.sign = 0;
X z2.sign = 0;
X if (ziszero(z2)) { /* number raised to power 0 */
X if (ziszero(z1))
X math_error("Zero raised to zero power");
X *res = _one_;
X return;
X }
X if (zisleone(z1)) { /* 0, 1, or -1 raised to a power */
X ans = _one_;
X ans.sign = (BOOL)sign;
X if (*z1.v == 0)
X ans = _zero_;
X *res = ans;
X return;
X }
X if (zisbig(z2))
X math_error("Raising to very large power");
X power = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
X if (zistwo(z1)) { /* two raised to a power */
X zbitvalue((long) power, res);
X return;
X }
X /*
X * See if this is a power of ten
X */
X if (zistiny(z1) && (*z1.v == 10)) {
X ztenpow((long) power, res);
X res->sign = (BOOL)sign;
X return;
X }
X /*
X * Handle low powers specially
X */
X if (power <= 4) {
X switch ((int) power) {
X case 1:
X ans.len = z1.len;
X ans.v = alloc(ans.len);
X zcopyval(z1, ans);
X ans.sign = (BOOL)sign;
X *res = ans;
X return;
X case 2:
X zsquare(z1, res);
X return;
X case 3:
X zsquare(z1, &temp);
X zmul(z1, temp, res);
X zfree(temp);
X res->sign = (BOOL)sign;
X return;
X case 4:
X zsquare(z1, &temp);
X zsquare(temp, res);
X zfree(temp);
X return;
X }
X }
X /*
X * Shift out all powers of twos so the multiplies are smaller.
X * We will shift back the right amount when done.
X */
X twos = 0;
X if (ziseven(z1)) {
X twos = zlowbit(z1);
X ans.v = alloc(z1.len);
X ans.len = z1.len;
X zcopyval(z1, ans);
X zshiftr(ans, twos);
X ztrim(&ans);
X z1 = ans;
X twos *= power;
X }
X /*
X * Compute the power by squaring and multiplying.
X * This uses the left to right method of power raising.
X */
X bit = TOPFULL;
X while ((bit & power) == 0)
X bit >>= 1L;
X bit >>= 1L;
X zsquare(z1, &ans);
X if (bit & power) {
X zmul(ans, z1, &temp);
X zfree(ans);
X ans = temp;
X }
X bit >>= 1L;
X while (bit) {
X zsquare(ans, &temp);
X zfree(ans);
X ans = temp;
X if (bit & power) {
X zmul(ans, z1, &temp);
X zfree(ans);
X ans = temp;
X }
X bit >>= 1L;
X }
X /*
X * Scale back up by proper power of two
X */
X if (twos) {
X zshift(ans, twos, &temp);
X zfree(ans);
X ans = temp;
X zfree(z1);
X }
X ans.sign = (BOOL)sign;
X *res = ans;
X}
X
X
X/*
X * Compute ten to the specified power
X * This saves some work since the squares of ten are saved.
X */
Xvoid
Xztenpow(power, res)
X long power;
X ZVALUE *res;
X{
X long i;
X ZVALUE ans;
X ZVALUE temp;
X
X if (power <= 0) {
X *res = _one_;
X return;
X }
X ans = _one_;
X _tenpowers_[0] = _ten_;
X for (i = 0; power; i++) {
X if (_tenpowers_[i].len == 0)
X zsquare(_tenpowers_[i-1], &_tenpowers_[i]);
X if (power & 0x1) {
X zmul(ans, _tenpowers_[i], &temp);
X zfree(ans);
X ans = temp;
X }
X power /= 2;
X }
X *res = ans;
X}
X
X
X/*
X * Calculate modular inverse suppressing unnecessary divisions.
X * This is based on the Euclidian algorithm for large numbers.
X * (Algorithm X from Knuth Vol 2, section 4.5.2. and exercise 17)
X * Returns TRUE if there is no solution because the numbers
X * are not relatively prime.
X */
XBOOL
Xzmodinv(u, v, res)
X ZVALUE u, v;
X ZVALUE *res;
X{
X FULL q1, q2, ui3, vi3, uh, vh, A, B, C, D, T;
X ZVALUE u2, u3, v2, v3, qz, tmp1, tmp2, tmp3;
X
X if (zisneg(u) || zisneg(v) || (zrel(u, v) >= 0))
X zmod(u, v, &v3);
X else
X zcopy(u, &v3);
X zcopy(v, &u3);
X u2 = _zero_;
X v2 = _one_;
X
X /*
X * Loop here while the size of the numbers remain above
X * the size of a FULL. Throughout this loop u3 >= v3.
X */
X while ((u3.len > 1) && !ziszero(v3)) {
X uh = (((FULL) u3.v[u3.len - 1]) << BASEB) + u3.v[u3.len - 2];
X vh = 0;
X if ((v3.len + 1) >= u3.len)
X vh = v3.v[v3.len - 1];
X if (v3.len == u3.len)
X vh = (vh << BASEB) + v3.v[v3.len - 2];
X A = 1;
X B = 0;
X C = 0;
X D = 1;
X
X /*
X * Calculate successive quotients of the continued fraction
X * expansion using only single precision arithmetic until
X * greater precision is required.
X */
X while ((vh + C) && (vh + D)) {
X q1 = (uh + A) / (vh + C);
X q2 = (uh + B) / (vh + D);
X if (q1 != q2)
X break;
X T = A - q1 * C;
X A = C;
X C = T;
X T = B - q1 * D;
X B = D;
X D = T;
X T = uh - q1 * vh;
X uh = vh;
X vh = T;
X }
X
X /*
X * If B is zero, then we made no progress because
X * the calculation requires a very large quotient.
X * So we must do this step of the calculation in
X * full precision
X */
X if (B == 0) {
X zquo(u3, v3, &qz);
X zmul(qz, v2, &tmp1);
X zsub(u2, tmp1, &tmp2);
X zfree(tmp1);
X zfree(u2);
X u2 = v2;
X v2 = tmp2;
X zmul(qz, v3, &tmp1);
X zsub(u3, tmp1, &tmp2);
X zfree(tmp1);
X zfree(u3);
X u3 = v3;
X v3 = tmp2;
X zfree(qz);
X continue;
X }
X /*
X * Apply the calculated A,B,C,D numbers to the current
X * values to update them as if the full precision
X * calculations had been carried out.
X */
X zmuli(u2, (long) A, &tmp1);
X zmuli(v2, (long) B, &tmp2);
X zadd(tmp1, tmp2, &tmp3);
X zfree(tmp1);
X zfree(tmp2);
X zmuli(u2, (long) C, &tmp1);
X zmuli(v2, (long) D, &tmp2);
X zfree(u2);
X zfree(v2);
X u2 = tmp3;
X zadd(tmp1, tmp2, &v2);
X zfree(tmp1);
X zfree(tmp2);
X zmuli(u3, (long) A, &tmp1);
X zmuli(v3, (long) B, &tmp2);
X zadd(tmp1, tmp2, &tmp3);
X zfree(tmp1);
X zfree(tmp2);
X zmuli(u3, (long) C, &tmp1);
X zmuli(v3, (long) D, &tmp2);
X zfree(u3);
X zfree(v3);
X u3 = tmp3;
X zadd(tmp1, tmp2, &v3);
X zfree(tmp1);
X zfree(tmp2);
X }
X
X /*
X * Here when the remaining numbers become single precision in size.
X * Finish the procedure using single precision calculations.
X */
X if (ziszero(v3) && !zisone(u3)) {
X zfree(u3);
X zfree(v3);
X zfree(u2);
X zfree(v2);
X return TRUE;
X }
X ui3 = (zistiny(u3) ? z1tol(u3) : z2tol(u3));
X vi3 = (zistiny(v3) ? z1tol(v3) : z2tol(v3));
X zfree(u3);
X zfree(v3);
X while (vi3) {
X q1 = ui3 / vi3;
X zmuli(v2, (long) q1, &tmp1);
X zsub(u2, tmp1, &tmp2);
X zfree(tmp1);
X zfree(u2);
X u2 = v2;
X v2 = tmp2;
X q2 = ui3 - q1 * vi3;
X ui3 = vi3;
X vi3 = q2;
X }
X zfree(v2);
X if (ui3 != 1) {
X zfree(u2);
X return TRUE;
X }
X if (zisneg(u2)) {
X zadd(v, u2, res);
X zfree(u2);
X return FALSE;
X }
X *res = u2;
X return FALSE;
X}
X
X
X#if 0
X/*
X * Approximate the quotient of two integers by another set of smaller
X * integers. This uses continued fractions to determine the smaller set.
X */
Xvoid
Xzapprox(z1, z2, res1, res2)
X ZVALUE z1, z2, *res1, *res2;
X{
X int sign;
X ZVALUE u1, v1, u3, v3, q, t1, t2, t3;
X
X sign = ((z1.sign != 0) ^ (z2.sign != 0));
X z1.sign = 0;
X z2.sign = 0;
X v3 = z2;
X u3 = z1;
X u1 = _one_;
X v1 = _zero_;
X while (!ziszero(v3)) {
X zdiv(u3, v3, &q, &t1);
X zmul(v1, q, &t2);
X zsub(u1, t2, &t3);
X zfree(q);
X zfree(t2);
X zfree(u1);
X if ((u3.v != z1.v) && (u3.v != z2.v))
X zfree(u3);
X u1 = v1;
X u3 = v3;
X v1 = t3;
X v3 = t1;
X }
X if (!zisunit(u3))
X math_error("Non-relativly prime numbers for approx");
X if ((u3.v != z1.v) && (u3.v != z2.v))
X zfree(u3);
X if ((v3.v != z1.v) && (v3.v != z2.v))
X zfree(v3);
X zfree(v1);
X zmul(u1, z1, &t1);
X zsub(t1, _one_, &t2);
X zfree(t1);
X zquo(t2, z2, &t1);
X zfree(t2);
X u1.sign = (BOOL)sign;
X t1.sign = 0;
X *res1 = t1;
X *res2 = u1;
X}
X#endif
X
X
X/*
X * Binary gcd algorithm
X * This algorithm taken from Knuth
X */
Xvoid
Xzgcd(z1, z2, res)
X ZVALUE z1, z2, *res;
X{
X ZVALUE u, v, t;
X register long j, k, olen, mask;
X register HALF h;
X HALF *oldv1, *oldv2;
X
X z1.sign = 0;
X z2.sign = 0;
X if (ziszero(z1)) {
X zcopy(z2, res);
X return;
X }
X if (ziszero(z2)) {
X zcopy(z1, res);
X return;
X }
X if (zisunit(z1) || zisunit(z2)) {
X *res = _one_;
X return;
X }
X /*
X * First see if one number is very much larger than the other.
X * If so, then divide as necessary to get the numbers near each other.
X */
X oldv1 = z1.v;
X oldv2 = z2.v;
X if (z1.len < z2.len) {
X t = z1;
X z1 = z2;
X z2 = t;
X }
X while ((z1.len > (z2.len + 5)) && !ziszero(z2)) {
X zmod(z1, z2, &t);
X if ((z1.v != oldv1) && (z1.v != oldv2))
X zfree(z1);
X z1 = z2;
X z2 = t;
X }
X /*
X * Ok, now do the binary method proper
X */
X u.len = z1.len;
X v.len = z2.len;
X u.sign = 0;
X v.sign = 0;
X if (!ztest(z1)) {
X v.v = alloc(v.len);
X zcopyval(z2, v);
X *res = v;
X goto done;
X }
X if (!ztest(z2)) {
X u.v = alloc(u.len);
X zcopyval(z1, u);
X *res = u;
X goto done;
X }
X u.v = alloc(u.len);
X v.v = alloc(v.len);
X zcopyval(z1, u);
X zcopyval(z2, v);
X k = 0;
X while (u.v[k] == 0 && v.v[k] == 0)
X ++k;
X mask = 01;
X h = u.v[k] | v.v[k];
X k *= BASEB;
X while (!(h & mask)) {
X mask <<= 1;
X k++;
X }
X zshiftr(u, k);
X zshiftr(v, k);
X ztrim(&u);
X ztrim(&v);
X if (zisodd(u)) {
X t.v = alloc(v.len);
X t.len = v.len;
X zcopyval(v, t);
X t.sign = !v.sign;
X } else {
X t.v = alloc(u.len);
X t.len = u.len;
X zcopyval(u, t);
X t.sign = u.sign;
X }
X while (ztest(t)) {
X j = 0;
X while (!t.v[j])
X ++j;
X mask = 01;
X h = t.v[j];
X j *= BASEB;
X while (!(h & mask)) {
X mask <<= 1;
X j++;
X }
X zshiftr(t, j);
X ztrim(&t);
X if (ztest(t) > 0) {
X zfree(u);
X u = t;
X } else {
X zfree(v);
X v = t;
X v.sign = !t.sign;
X }
X zsub(u, v, &t);
X }
X zfree(t);
X zfree(v);
X if (k) {
X olen = u.len;
X u.len += k / BASEB + 1;
X u.v = (HALF *) realloc(u.v, u.len * sizeof(HALF));
X while (olen != u.len)
X u.v[olen++] = 0;
X zshiftl(u, k);
X }
X ztrim(&u);
X *res = u;
X
Xdone:
X if ((z1.v != oldv1) && (z1.v != oldv2))
X zfree(z1);
X if ((z2.v != oldv1) && (z2.v != oldv2))
X zfree(z2);
X}
X
X
X/*
X * Compute the lcm of two integers (least common multiple).
X * This is done using the formula: gcd(a,b) * lcm(a,b) = a * b.
X */
Xvoid
Xzlcm(z1, z2, res)
X ZVALUE z1, z2, *res;
X{
X ZVALUE temp1, temp2;
X
X zgcd(z1, z2, &temp1);
X zquo(z1, temp1, &temp2);
X zfree(temp1);
X zmul(temp2, z2, res);
X zfree(temp2);
X}
X
X
X/*
X * Return whether or not two numbers are relatively prime to each other.
X */
XBOOL
Xzrelprime(z1, z2)
X ZVALUE z1, z2; /* numbers to be tested */
X{
X FULL rem1, rem2; /* remainders */
X ZVALUE rem;
X BOOL result;
X
X z1.sign = 0;
X z2.sign = 0;
X if (ziseven(z1) && ziseven(z2)) /* false if both even */
X return FALSE;
X if (zisunit(z1) || zisunit(z2)) /* true if either is a unit */
X return TRUE;
X if (ziszero(z1) || ziszero(z2)) /* false if either is zero */
X return FALSE;
X if (zistwo(z1) || zistwo(z2)) /* true if either is two */
X return TRUE;
X /*
X * Try reducing each number by the product of the first few odd primes
X * to see if any of them are a common factor.
X */
X rem1 = zmodi(z1, 3L * 5 * 7 * 11 * 13);
X rem2 = zmodi(z2, 3L * 5 * 7 * 11 * 13);
X if (((rem1 % 3) == 0) && ((rem2 % 3) == 0))
X return FALSE;
X if (((rem1 % 5) == 0) && ((rem2 % 5) == 0))
X return FALSE;
X if (((rem1 % 7) == 0) && ((rem2 % 7) == 0))
X return FALSE;
X if (((rem1 % 11) == 0) && ((rem2 % 11) == 0))
X return FALSE;
X if (((rem1 % 13) == 0) && ((rem2 % 13) == 0))
X return FALSE;
X /*
X * Try a new batch of primes now
X */
X rem1 = zmodi(z1, 17L * 19 * 23);
X rem2 = zmodi(z2, 17L * 19 * 23);
X if (((rem1 % 17) == 0) && ((rem2 % 17) == 0))
X return FALSE;
X if (((rem1 % 19) == 0) && ((rem2 % 19) == 0))
X return FALSE;
X if (((rem1 % 23) == 0) && ((rem2 % 23) == 0))
X return FALSE;
X /*
X * Yuk, we must actually compute the gcd to know the answer
X */
X zgcd(z1, z2, &rem);
X result = zisunit(rem);
X zfree(rem);
X return result;
X}
X
X
X/*
X * Compute the log of one number base another, to the closest integer.
X * This is the largest integer which when the second number is raised to it,
X * the resulting value is less than or equal to the first number.
X * Example: zlog(123456, 10) = 5.
X */
Xlong
Xzlog(z1, z2)
X ZVALUE z1, z2;
X{
X register ZVALUE *zp; /* current square */
X long power; /* current power */
X long worth; /* worth of current square */
X ZVALUE val; /* current value of power */
X ZVALUE temp; /* temporary */
X ZVALUE squares[32]; /* table of squares of base */
X
X /*
X * Make sure that the numbers are nonzero and the base is greater than one.
X */
X if (zisneg(z1) || ziszero(z1) || zisneg(z2) || zisleone(z2))
X math_error("Bad arguments for log");
X /*
X * Reject trivial cases.
X */
X if (z1.len < z2.len)
X return 0;
X if ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))
X return 0;
X power = zrel(z1, z2);
X if (power <= 0)
X return (power + 1);
X /*
X * Handle any power of two special.
X */
X if (zisonebit(z2))
X return (zhighbit(z1) / zlowbit(z2));
X /*
X * Handle base 10 special
X */
X if ((z2.len == 1) && (*z2.v == 10))
X return zlog10(z1);
X /*
X * Now loop by squaring the base each time, and see whether or
X * not each successive square is still smaller than the number.
X */
X worth = 1;
X zp = &squares[0];
X *zp = z2;
X while (((zp->len * 2) - 1) <= z1.len) { /* while square not too large */
X zsquare(*zp, zp + 1);
X zp++;
X worth *= 2;
X }
X /*
X * Now back down the squares, and multiply them together to see
X * exactly how many times the base can be raised by.
X */
X val = _one_;
X power = 0;
X for (; zp >= squares; zp--, worth /= 2) {
X if ((val.len + zp->len - 1) <= z1.len) {
X zmul(val, *zp, &temp);
X if (zrel(z1, temp) >= 0) {
X zfree(val);
X val = temp;
X power += worth;
X } else
X zfree(temp);
X }
X if (zp != squares)
X zfree(*zp);
X }
X return power;
X}
X
X
X/*
X * Return the integral log base 10 of a number.
X */
Xlong
Xzlog10(z)
X ZVALUE z;
X{
X register ZVALUE *zp; /* current square */
X long power; /* current power */
X long worth; /* worth of current square */
X ZVALUE val; /* current value of power */
X ZVALUE temp; /* temporary */
X
X if (!zispos(z))
X math_error("Non-positive number for log10");
X /*
X * Loop by squaring the base each time, and see whether or
X * not each successive square is still smaller than the number.
X */
X worth = 1;
X zp = &_tenpowers_[0];
X *zp = _ten_;
X while (((zp->len * 2) - 1) <= z.len) { /* while square not too large */
X if (zp[1].len == 0)
X zsquare(*zp, zp + 1);
X zp++;
X worth *= 2;
X }
X /*
X * Now back down the squares, and multiply them together to see
X * exactly how many times the base can be raised by.
X */
X val = _one_;
X power = 0;
X for (; zp >= _tenpowers_; zp--, worth /= 2) {
X if ((val.len + zp->len - 1) <= z.len) {
X zmul(val, *zp, &temp);
X if (zrel(z, temp) >= 0) {
X zfree(val);
X val = temp;
X power += worth;
X } else
X zfree(temp);
X }
X }
X return power;
X}
X
X
X/*
X * Return the number of times that one number will divide another.
X * This works similarly to zlog, except that divisions must be exact.
X * For example, zdivcount(540, 3) = 3, since 3^3 divides 540 but 3^4 won't.
X */
Xlong
Xzdivcount(z1, z2)
X ZVALUE z1, z2;
X{
X long count; /* number of factors removed */
X ZVALUE tmp; /* ignored return value */
X
X count = zfacrem(z1, z2, &tmp);
X zfree(tmp);
X return count;
X}
X
X
X/*
X * Remove all occurances of the specified factor from a number.
X * Also returns the number of factors removed as a function return value.
X * Example: zfacrem(540, 3, &x) returns 3 and sets x to 20.
X */
Xlong
Xzfacrem(z1, z2, rem)
X ZVALUE z1, z2, *rem;
X{
X register ZVALUE *zp; /* current square */
X long count; /* total count of divisions */
X long worth; /* worth of current square */
X ZVALUE temp1, temp2, temp3; /* temporaries */
X ZVALUE squares[32]; /* table of squares of factor */
X
X z1.sign = 0;
X z2.sign = 0;
X /*
X * Make sure factor isn't 0 or 1.
X */
X if (zisleone(z2))
X math_error("Bad argument for facrem");
X /*
X * Reject trivial cases.
X */
X if ((z1.len < z2.len) || (zisodd(z1) && ziseven(z2)) ||
X ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))) {
X rem->v = alloc(z1.len);
X rem->len = z1.len;
X rem->sign = 0;
X zcopyval(z1, *rem);
X return 0;
X }
X /*
X * Handle any power of two special.
X */
X if (zisonebit(z2)) {
X count = zlowbit(z1) / zlowbit(z2);
X rem->v = alloc(z1.len);
X rem->len = z1.len;
X rem->sign = 0;
X zcopyval(z1, *rem);
X zshiftr(*rem, count);
X ztrim(rem);
X return count;
X }
X /*
X * See if the factor goes in even once.
X */
X zdiv(z1, z2, &temp1, &temp2);
X if (!ziszero(temp2)) {
X zfree(temp1);
X zfree(temp2);
X rem->v = alloc(z1.len);
X rem->len = z1.len;
X rem->sign = 0;
X zcopyval(z1, *rem);
X return 0;
X }
X zfree(temp2);
X z1 = temp1;
X /*
X * Now loop by squaring the factor each time, and see whether
X * or not each successive square will still divide the number.
X */
X count = 1;
X worth = 1;
X zp = &squares[0];
X *zp = z2;
X while (((zp->len * 2) - 1) <= z1.len) { /* while square not too large */
X zsquare(*zp, &temp1);
X zdiv(z1, temp1, &temp2, &temp3);
X if (!ziszero(temp3)) {
X zfree(temp1);
X zfree(temp2);
X zfree(temp3);
X break;
X }
X zfree(temp3);
X zfree(z1);
X z1 = temp2;
X *++zp = temp1;
X worth *= 2;
X count += worth;
X }
X /*
X * Now back down the list of squares, and see if the lower powers
X * will divide any more times.
X */
X for (; zp >= squares; zp--, worth /= 2) {
X if (zp->len <= z1.len) {
X zdiv(z1, *zp, &temp1, &temp2);
X if (ziszero(temp2)) {
X temp3 = z1;
X z1 = temp1;
X temp1 = temp3;
X count += worth;
X }
X zfree(temp1);
X zfree(temp2);
X }
X if (zp != squares)
X zfree(*zp);
X }
X *rem = z1;
X return count;
X}
X
X
X/*
X * Keep dividing a number by the gcd of it with another number until the
X * result is relatively prime to the second number.
X */
Xvoid
Xzgcdrem(z1, z2, res)
X ZVALUE z1, z2, *res;
X{
X ZVALUE tmp1, tmp2;
X
X /*
X * Begin by taking the gcd for the first time.
X * If the number is already relatively prime, then we are done.
X */
X zgcd(z1, z2, &tmp1);
X if (zisunit(tmp1) || ziszero(tmp1)) {
X res->len = z1.len;
X res->v = alloc(z1.len);
X res->sign = z1.sign;
X zcopyval(z1, *res);
X return;
X }
X zquo(z1, tmp1, &tmp2);
X z1 = tmp2;
X z2 = tmp1;
X /*
X * Now keep alternately taking the gcd and removing factors until
X * the gcd becomes one.
X */
X while (!zisunit(z2)) {
X (void) zfacrem(z1, z2, &tmp1);
X zfree(z1);
X z1 = tmp1;
X zgcd(z1, z2, &tmp1);
X zfree(z2);
X z2 = tmp1;
X }
X *res = z1;
X}
X
X
X/*
X * Find the lowest prime factor of a number if one can be found.
X * Search is conducted for the first count primes.
X * Returns 1 if no factor was found.
X */
Xlong
Xzlowfactor(z, count)
X ZVALUE z;
X long count;
X{
X FULL p, d;
X ZVALUE div, tmp;
X HALF divval[2];
X
X if ((--count < 0) || ziszero(z))
X return 1;
X if (ziseven(z))
X return 2;
X div.sign = 0;
X div.v = divval;
X for (p = 3; (count > 0); p += 2) {
X for (d = 3; (d * d) <= p; d += 2)
X if ((p % d) == 0)
X goto next;
X divval[0] = (p & BASE1);
X divval[1] = (p / BASE);
X div.len = 1 + (p >= BASE);
X zmod(z, div, &tmp);
X if (ziszero(tmp))
X return p;
X zfree(tmp);
X count--;
Xnext:;
X }
X return 1;
X}
X
X
X/*
X * Return the number of digits (base 10) in a number, ignoring the sign.
X */
Xlong
Xzdigits(z1)
X ZVALUE z1;
X{
X long count, val;
X
X z1.sign = 0;
X if (zistiny(z1)) { /* do small numbers ourself */
X count = 1;
X val = 10;
X while (*z1.v >= (HALF)val) {
X count++;
X val *= 10;
X }
X return count;
X }
X return (zlog10(z1) + 1);
X}
X
X
X/*
X * Return the single digit at the specified decimal place of a number,
X * where 0 means the rightmost digit. Example: zdigit(1234, 1) = 3.
X */
XFLAG
Xzdigit(z1, n)
X ZVALUE z1;
X long n;
X{
X ZVALUE tmp1, tmp2;
X FLAG res;
X
X z1.sign = 0;
X if (ziszero(z1) || (n < 0) || (n / BASEDIG >= z1.len))
X return 0;
X if (n == 0)
X return zmodi(z1, 10L);
X if (n == 1)
X return zmodi(z1, 100L) / 10;
X if (n == 2)
X return zmodi(z1, 1000L) / 100;
X if (n == 3)
X return zmodi(z1, 10000L) / 1000;
X ztenpow(n, &tmp1);
X zquo(z1, tmp1, &tmp2);
X res = zmodi(tmp2, 10L);
X zfree(tmp1);
X zfree(tmp2);
X return res;
X}
X
X
X/*
X * Find the square root of a number. This is the greatest integer whose
X * square is less than or equal to the number. Returns TRUE if the
X * square root is exact.
X */
XBOOL
Xzsqrt(z1, dest)
X ZVALUE z1, *dest;
X{
X ZVALUE try, quo, rem, old, temp;
X FULL iquo, val;
X long i,j;
X
X if (z1.sign)
X math_error("Square root of negative number");
X if (ziszero(z1)) {
X *dest = _zero_;
X return TRUE;
X }
X if ((*z1.v < 4) && zistiny(z1)) {
X *dest = _one_;
X return (*z1.v == 1);
X }
X /*
X * Pick the square root of the leading one or two digits as a first guess.
X */
X val = z1.v[z1.len-1];
X if ((z1.len & 0x1) == 0)
X val = (val * BASE) + z1.v[z1.len-2];
X
X /*
X * Find the largest power of 2 that when squared
X * is <= val > 0. Avoid multiply overflow by doing
X * a careful check at the BASE boundary.
X */
X j = 1<<(BASEB+BASEB-2);
X if (val > j) {
X iquo = BASE;
X } else {
X i = BASEB-1;
X while (j > val) {
X --i;
X j >>= 2;
X }
X iquo = bitmask[i];
X }
X
X for (i = 8; i > 0; i--)
X iquo = (iquo + (val / iquo)) / 2;
X if (iquo > BASE1)
X iquo = BASE1;
X /*
X * Allocate the numbers to use for the main loop.
X * The size and high bits of the final result are correctly set here.
X * Notice that the remainder of the test value is rubbish, but this
X * is unimportant.
X */
X try.sign = 0;
X try.len = (z1.len + 1) / 2;
X try.v = alloc(try.len);
X try.v[try.len-1] = (HALF)iquo;
X old.sign = 0;
X old.v = alloc(try.len);
X old.len = 1;
X *old.v = 0;
X /*
X * Main divide and average loop
X */
X for (;;) {
X zdiv(z1, try, &quo, &rem);
X i = zrel(try, quo);
X if ((i == 0) && ziszero(rem)) { /* exact square root */
X zfree(rem);
X zfree(quo);
X zfree(old);
X *dest = try;
X return TRUE;
X }
X zfree(rem);
X if (i <= 0) {
X /*
X * Current try is less than or equal to the square root since
X * it is less than the quotient. If the quotient is equal to
X * the try, we are all done. Also, if the try is equal to the
X * old try value, we are done since no improvement occurred.
X * If not, save the improved value and loop some more.
X */
X if ((i == 0) || (zcmp(old, try) == 0)) {
X zfree(quo);
X zfree(old);
X *dest = try;
X return FALSE;
X }
X old.len = try.len;
X zcopyval(try, old);
X }
X /* average current try and quotent for the new try */
X zadd(try, quo, &temp);
X zfree(quo);
X zfree(try);
X try = temp;
X zshiftr(try, 1L);
X zquicktrim(try);
X }
X}
X
X
X/*
X * Take an arbitrary root of a number (to the greatest integer).
X * This uses the following iteration to get the Kth root of N:
X * x = ((K-1) * x + N / x^(K-1)) / K
X */
Xvoid
Xzroot(z1, z2, dest)
X ZVALUE z1, z2, *dest;
X{
X ZVALUE try, quo, old, temp, temp2;
X ZVALUE k1; /* holds k - 1 */
X int sign;
X long i, k, highbit;
X SIUNION sival;
X
X sign = z1.sign;
X if (sign && ziseven(z2))
X math_error("Even root of negative number");
X if (ziszero(z2) || zisneg(z2))
X math_error("Non-positive root");
X if (ziszero(z1)) { /* root of zero */
X *dest = _zero_;
X return;
X }
X if (zisunit(z2)) { /* first root */
X zcopy(z1, dest);
X return;
X }
X if (zisbig(z2)) { /* humongous root */
X *dest = _one_;
X dest->sign = (HALF)sign;
X return;
X }
X k = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
X highbit = zhighbit(z1);
X if (highbit < k) { /* too high a root */
X *dest = _one_;
X dest->sign = (HALF)sign;
X return;
X }
X sival.ivalue = k - 1;
X k1.v = &sival.silow;
X k1.len = 1 + (sival.sihigh != 0);
X k1.sign = 0;
X z1.sign = 0;
X /*
X * Allocate the numbers to use for the main loop.
X * The size and high bits of the final result are correctly set here.
X * Notice that the remainder of the test value is rubbish, but this
SHAR_EOF
echo "End of part 11"
echo "File calc2.9.0/zfunc.c is continued in part 12"
echo "12" > s2_seq_.tmp
exit 0
