BSD 4_1_snap development

Work on file usr/src/lib/libm/atan.c
Work on file usr/src/lib/libm/floor.c
Work on file usr/src/lib/libm/hypot.c
Work on file usr/src/lib/libm/erf.c
Work on file usr/src/lib/libm/exp.c
Work on file usr/src/lib/libm/sinh.c
Work on file usr/src/lib/libm/tan.c
Work on file usr/src/lib/libmp/madd.c
Work on file usr/src/lib/libm/pow.c
Work on file usr/src/lib/libm/jn.c
Work on file usr/src/lib/libm/sin.c
Work on file usr/src/lib/libm/j0.c
Work on file usr/src/lib/libnm/pow.c
Work on file usr/src/lib/libmp/gcd.c
Work on file usr/src/lib/libm/j1.c
Work on file usr/src/lib/libmp/mout.c
Work on file usr/src/lib/libmp/mdiv.c
Work on file usr/src/lib/libm/sqrt.c
Work on file usr/src/lib/libm/tanh.c
Work on file usr/src/lib/libm/log.c
Work on file usr/src/lib/libmp/mult.c
Work on file usr/src/lib/libmp/msqrt.c
Work on file usr/src/lib/libmp/util.c
Work on file usr/src/lib/libmp/pow.c

Synthesized-from: CSRG/cd1/4.1.snap

diff --git a/usr/src/lib/libm/atan.c b/usr/src/lib/libm/atan.c
new file mode 100644 (file)
index 0000000..3c8d07e
--- /dev/null
+++ b/usr/src/lib/libm/atan.c
@@ -0,0 +1,111 @@
+
+/*
+       floating-point arctangent
+
+       atan returns the value of the arctangent of its
+       argument in the range [-pi/2,pi/2].
+
+       atan2 returns the arctangent of arg1/arg2
+       in the range [-pi,pi].
+
+       there are no error returns.
+
+       coefficients are #5077 from Hart & Cheney. (19.56D)
+*/
+
+
+double static sq2p1     =2.414213562373095048802e0;
+static double sq2m1     = .414213562373095048802e0;
+static double pio2      =1.570796326794896619231e0;
+static double pio4      = .785398163397448309615e0;
+static double p4        = .161536412982230228262e2;
+static double p3        = .26842548195503973794141e3;
+static double p2        = .11530293515404850115428136e4;
+static double p1        = .178040631643319697105464587e4;
+static double p0        = .89678597403663861959987488e3;
+static double q4        = .5895697050844462222791e2;
+static double q3        = .536265374031215315104235e3;
+static double q2        = .16667838148816337184521798e4;
+static double q1        = .207933497444540981287275926e4;
+static double q0        = .89678597403663861962481162e3;
+
+
+/*
+       atan makes its argument positive and
+       calls the inner routine satan.
+*/
+
+double
+atan(arg)
+double arg;
+{
+       double satan();
+
+       if(arg>0)
+               return(satan(arg));
+       else
+               return(-satan(-arg));
+}
+
+
+/*
+       atan2 discovers what quadrant the angle
+       is in and calls atan.
+*/
+
+double
+atan2(arg1,arg2)
+double arg1,arg2;
+{
+       double satan();
+
+       if((arg1+arg2)==arg1)
+               if(arg1 >= 0.) return(pio2);
+               else return(-pio2);
+       else if(arg2 <0.)
+               if(arg1 >= 0.)
+                       return(pio2+pio2 - satan(-arg1/arg2));
+               else
+                       return(-pio2-pio2 + satan(arg1/arg2));
+       else if(arg1>0)
+               return(satan(arg1/arg2));
+       else
+               return(-satan(-arg1/arg2));
+}
+
+/*
+       satan reduces its argument (known to be positive)
+       to the range [0,0.414...] and calls xatan.
+*/
+
+static double
+satan(arg)
+double arg;
+{
+       double  xatan();
+
+       if(arg < sq2m1)
+               return(xatan(arg));
+       else if(arg > sq2p1)
+               return(pio2 - xatan(1.0/arg));
+       else
+               return(pio4 + xatan((arg-1.0)/(arg+1.0)));
+}
+
+/*
+       xatan evaluates a series valid in the
+       range [-0.414...,+0.414...].
+*/
+
+static double
+xatan(arg)
+double arg;
+{
+       double argsq;
+       double value;
+
+       argsq = arg*arg;
+       value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0);
+       value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0);
+       return(value*arg);
+}

diff --git a/usr/src/lib/libm/erf.c b/usr/src/lib/libm/erf.c
new file mode 100644 (file)
index 0000000..b724b90
--- /dev/null
+++ b/usr/src/lib/libm/erf.c
@@ -0,0 +1,116 @@
+/*
+       C program for floating point error function
+
+       erf(x) returns the error function of its argument
+       erfc(x) returns 1.0-erf(x)
+
+       erf(x) is defined by
+       ${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$
+
+       the entry for erfc is provided because of the
+       extreme loss of relative accuracy if erf(x) is
+       called for large x and the result subtracted
+       from 1. (e.g. for x= 10, 12 places are lost).
+
+       There are no error returns.
+
+       Calls exp.
+
+       Coefficients for large x are #5667 from Hart & Cheney (18.72D).
+*/
+
+#define M 7
+#define N 9
+int errno;
+static double torp = 1.1283791670955125738961589031;
+static double p1[] = {
+       0.804373630960840172832162e5,
+       0.740407142710151470082064e4,
+       0.301782788536507577809226e4,
+       0.380140318123903008244444e2,
+       0.143383842191748205576712e2,
+       -.288805137207594084924010e0,
+       0.007547728033418631287834e0,
+};
+static double q1[]  = {
+       0.804373630960840172826266e5,
+       0.342165257924628539769006e5,
+       0.637960017324428279487120e4,
+       0.658070155459240506326937e3,
+       0.380190713951939403753468e2,
+       0.100000000000000000000000e1,
+       0.0,
+};
+static double p2[]  = {
+       0.18263348842295112592168999e4,
+       0.28980293292167655611275846e4,
+       0.2320439590251635247384768711e4,
+       0.1143262070703886173606073338e4,
+       0.3685196154710010637133875746e3,
+       0.7708161730368428609781633646e2,
+       0.9675807882987265400604202961e1,
+       0.5641877825507397413087057563e0,
+       0.0,
+};
+static double q2[]  = {
+       0.18263348842295112595576438e4,
+       0.495882756472114071495438422e4,
+       0.60895424232724435504633068e4,
+       0.4429612803883682726711528526e4,
+       0.2094384367789539593790281779e4,
+       0.6617361207107653469211984771e3,
+       0.1371255960500622202878443578e3,
+       0.1714980943627607849376131193e2,
+       1.0,
+};
+
+double
+erf(arg) double arg;{
+       double erfc();
+       int sign;
+       double argsq;
+       double d, n;
+       int i;
+
+       errno = 0;
+       sign = 1;
+       if(arg < 0.){
+               arg = -arg;
+               sign = -1;
+       }
+       if(arg < 0.5){
+               argsq = arg*arg;
+               for(n=0,d=0,i=M-1; i>=0; i--){
+                       n = n*argsq + p1[i];
+                       d = d*argsq + q1[i];
+               }
+               return(sign*torp*arg*n/d);
+       }
+       if(arg >= 10.)
+               return(sign*1.);
+       return(sign*(1. - erfc(arg)));
+}
+
+double
+erfc(arg) double arg;{
+       double erf();
+       double exp();
+       double n, d;
+       int i;
+
+       errno = 0;
+       if(arg < 0.)
+               return(2. - erfc(-arg));
+/*
+       if(arg < 0.5)
+               return(1. - erf(arg));
+*/
+       if(arg >= 10.)
+               return(0.);
+
+       for(n=0,d=0,i=N-1; i>=0; i--){
+               n = n*arg + p2[i];
+               d = d*arg + q2[i];
+       }
+       return(exp(-arg*arg)*n/d);
+}

diff --git a/usr/src/lib/libm/exp.c b/usr/src/lib/libm/exp.c
new file mode 100644 (file)
index 0000000..84d28f0
--- /dev/null
+++ b/usr/src/lib/libm/exp.c
@@ -0,0 +1,45 @@
+/*
+       exp returns the exponential function of its
+       floating-point argument.
+
+       The coefficients are #1069 from Hart and Cheney. (22.35D)
+*/
+
+#include <errno.h>
+#include <math.h>
+
+int    errno;
+static double  p0      = .2080384346694663001443843411e7;
+static double  p1      = .3028697169744036299076048876e5;
+static double  p2      = .6061485330061080841615584556e2;
+static double  q0      = .6002720360238832528230907598e7;
+static double  q1      = .3277251518082914423057964422e6;
+static double  q2      = .1749287689093076403844945335e4;
+static double  log2e   = 1.4426950408889634073599247;
+static double  sqrt2   = 1.4142135623730950488016887;
+static double  maxf    = 10000;
+
+double
+exp(arg)
+double arg;
+{
+       double fract;
+       double temp1, temp2, xsq;
+       int ent;
+
+       if(arg == 0.)
+               return(1.);
+       if(arg < -maxf)
+               return(0.);
+       if(arg > maxf) {
+               errno = ERANGE;
+               return(HUGE);
+       }
+       arg *= log2e;
+       ent = floor(arg);
+       fract = (arg-ent) - 0.5;
+       xsq = fract*fract;
+       temp1 = ((p2*xsq+p1)*xsq+p0)*fract;
+       temp2 = ((1.0*xsq+q2)*xsq+q1)*xsq + q0;
+       return(ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent));
+}

diff --git a/usr/src/lib/libm/floor.c b/usr/src/lib/libm/floor.c
new file mode 100644 (file)
index 0000000..6c4aa70
--- /dev/null
+++ b/usr/src/lib/libm/floor.c
@@ -0,0 +1,30 @@
+/*
+ * floor and ceil-- greatest integer <= arg
+ * (resp least >=)
+ */
+
+double modf();
+
+double
+floor(d)
+double d;
+{
+       double fract;
+
+       if (d<0.0) {
+               d = -d;
+               fract = modf(d, &d);
+               if (fract != 0.0)
+                       d += 1;
+               d = -d;
+       } else
+               modf(d, &d);
+       return(d);
+}
+
+double
+ceil(d)
+double d;
+{
+       return(-floor(-d));
+}

diff --git a/usr/src/lib/libm/hypot.c b/usr/src/lib/libm/hypot.c
new file mode 100644 (file)
index 0000000..73d009f
--- /dev/null
+++ b/usr/src/lib/libm/hypot.c
@@ -0,0 +1,41 @@
+/*
+ * sqrt(a^2 + b^2)
+ *     (but carefully)
+ */
+
+double sqrt();
+double
+hypot(a,b)
+double a,b;
+{
+       double t;
+       if(a<0) a = -a;
+       if(b<0) b = -b;
+       if(a > b) {
+               t = a;
+               a = b;
+               b = t;
+       }
+       if(b==0) return(0.);
+       a /= b;
+       /*
+        * pathological overflow possible
+        * in the next line.
+        */
+       return(b*sqrt(1. + a*a));
+}
+
+struct complex
+{
+       double  r;
+       double  i;
+};
+
+double
+cabs(arg)
+struct complex arg;
+{
+       double hypot();
+
+       return(hypot(arg.r, arg.i));
+}

diff --git a/usr/src/lib/libm/j0.c b/usr/src/lib/libm/j0.c
new file mode 100644 (file)
index 0000000..2b2eb54
--- /dev/null
+++ b/usr/src/lib/libm/j0.c
@@ -0,0 +1,187 @@
+/*
+       floating point Bessel's function
+       of the first and second kinds
+       of order zero
+
+       j0(x) returns the value of J0(x)
+       for all real values of x.
+
+       There are no error returns.
+       Calls sin, cos, sqrt.
+
+       There is a niggling bug in J0 which
+       causes errors up to 2e-16 for x in the
+       interval [-8,8].
+       The bug is caused by an inappropriate order
+       of summation of the series.  rhm will fix it
+       someday.
+
+       Coefficients are from Hart & Cheney.
+       #5849 (19.22D)
+       #6549 (19.25D)
+       #6949 (19.41D)
+
+       y0(x) returns the value of Y0(x)
+       for positive real values of x.
+       For x<=0, error number EDOM is set and a
+       large negative value is returned.
+
+       Calls sin, cos, sqrt, log, j0.
+
+       The values of Y0 have not been checked
+       to more than ten places.
+
+       Coefficients are from Hart & Cheney.
+       #6245 (18.78D)
+       #6549 (19.25D)
+       #6949 (19.41D)
+*/
+
+#include <math.h>
+#include <errno.h>
+
+int    errno;
+static double pzero, qzero;
+static double tpi      = .6366197723675813430755350535e0;
+static double pio4     = .7853981633974483096156608458e0;
+static double p1[] = {
+       0.4933787251794133561816813446e21,
+       -.1179157629107610536038440800e21,
+       0.6382059341072356562289432465e19,
+       -.1367620353088171386865416609e18,
+       0.1434354939140344111664316553e16,
+       -.8085222034853793871199468171e13,
+       0.2507158285536881945555156435e11,
+       -.4050412371833132706360663322e8,
+       0.2685786856980014981415848441e5,
+};
+static double q1[] = {
+       0.4933787251794133562113278438e21,
+       0.5428918384092285160200195092e19,
+       0.3024635616709462698627330784e17,
+       0.1127756739679798507056031594e15,
+       0.3123043114941213172572469442e12,
+       0.6699987672982239671814028660e9,
+       0.1114636098462985378182402543e7,
+       0.1363063652328970604442810507e4,
+       1.0
+};
+static double p2[] = {
+       0.5393485083869438325262122897e7,
+       0.1233238476817638145232406055e8,
+       0.8413041456550439208464315611e7,
+       0.2016135283049983642487182349e7,
+       0.1539826532623911470917825993e6,
+       0.2485271928957404011288128951e4,
+       0.0,
+};
+static double q2[] = {
+       0.5393485083869438325560444960e7,
+       0.1233831022786324960844856182e8,
+       0.8426449050629797331554404810e7,
+       0.2025066801570134013891035236e7,
+       0.1560017276940030940592769933e6,
+       0.2615700736920839685159081813e4,
+       1.0,
+};
+static double p3[] = {
+       -.3984617357595222463506790588e4,
+       -.1038141698748464093880530341e5,
+       -.8239066313485606568803548860e4,
+       -.2365956170779108192723612816e4,
+       -.2262630641933704113967255053e3,
+       -.4887199395841261531199129300e1,
+       0.0,
+};
+static double q3[] = {
+       0.2550155108860942382983170882e6,
+       0.6667454239319826986004038103e6,
+       0.5332913634216897168722255057e6,
+       0.1560213206679291652539287109e6,
+       0.1570489191515395519392882766e5,
+       0.4087714673983499223402830260e3,
+       1.0,
+};
+static double p4[] = {
+       -.2750286678629109583701933175e20,
+       0.6587473275719554925999402049e20,
+       -.5247065581112764941297350814e19,
+       0.1375624316399344078571335453e18,
+       -.1648605817185729473122082537e16,
+       0.1025520859686394284509167421e14,
+       -.3436371222979040378171030138e11,
+       0.5915213465686889654273830069e8,
+       -.4137035497933148554125235152e5,
+};
+static double q4[] = {
+       0.3726458838986165881989980e21,
+       0.4192417043410839973904769661e19,
+       0.2392883043499781857439356652e17,
+       0.9162038034075185262489147968e14,
+       0.2613065755041081249568482092e12,
+       0.5795122640700729537480087915e9,
+       0.1001702641288906265666651753e7,
+       0.1282452772478993804176329391e4,
+       1.0,
+};
+
+double
+j0(arg) double arg;{
+       double argsq, n, d;
+       double sin(), cos(), sqrt();
+       int i;
+
+       if(arg < 0.) arg = -arg;
+       if(arg > 8.){
+               asympt(arg);
+               n = arg - pio4;
+               return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)));
+       }
+       argsq = arg*arg;
+       for(n=0,d=0,i=8;i>=0;i--){
+               n = n*argsq + p1[i];
+               d = d*argsq + q1[i];
+       }
+       return(n/d);
+}
+
+double
+y0(arg) double arg;{
+       double argsq, n, d;
+       double sin(), cos(), sqrt(), log(), j0();
+       int i;
+
+       errno = 0;
+       if(arg <= 0.){
+               errno = EDOM;
+               return(-HUGE);
+       }
+       if(arg > 8.){
+               asympt(arg);
+               n = arg - pio4;
+               return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)));
+       }
+       argsq = arg*arg;
+       for(n=0,d=0,i=8;i>=0;i--){
+               n = n*argsq + p4[i];
+               d = d*argsq + q4[i];
+       }
+       return(n/d + tpi*j0(arg)*log(arg));
+}
+
+static
+asympt(arg) double arg;{
+       double zsq, n, d;
+       int i;
+       zsq = 64./(arg*arg);
+       for(n=0,d=0,i=6;i>=0;i--){
+               n = n*zsq + p2[i];
+               d = d*zsq + q2[i];
+       }
+       pzero = n/d;
+       for(n=0,d=0,i=6;i>=0;i--){
+               n = n*zsq + p3[i];
+               d = d*zsq + q3[i];
+       }
+       qzero = (8./arg)*(n/d);
+}

diff --git a/usr/src/lib/libm/j1.c b/usr/src/lib/libm/j1.c
new file mode 100644 (file)
index 0000000..13b397f
--- /dev/null
+++ b/usr/src/lib/libm/j1.c
@@ -0,0 +1,193 @@
+/*
+       floating point Bessel's function
+       of the first and second kinds
+       of order one
+
+       j1(x) returns the value of J1(x)
+       for all real values of x.
+
+       There are no error returns.
+       Calls sin, cos, sqrt.
+
+       There is a niggling bug in J1 which
+       causes errors up to 2e-16 for x in the
+       interval [-8,8].
+       The bug is caused by an inappropriate order
+       of summation of the series.  rhm will fix it
+       someday.
+
+       Coefficients are from Hart & Cheney.
+       #6050 (20.98D)
+       #6750 (19.19D)
+       #7150 (19.35D)
+
+       y1(x) returns the value of Y1(x)
+       for positive real values of x.
+       For x<=0, error number EDOM is set and a
+       large negative value is returned.
+
+       Calls sin, cos, sqrt, log, j1.
+
+       The values of Y1 have not been checked
+       to more than ten places.
+
+       Coefficients are from Hart & Cheney.
+       #6447 (22.18D)
+       #6750 (19.19D)
+       #7150 (19.35D)
+*/
+
+#include <math.h>
+#include <errno.h>
+
+int    errno;
+static double pzero, qzero;
+static double tpi      = .6366197723675813430755350535e0;
+static double pio4     = .7853981633974483096156608458e0;
+static double p1[] = {
+       0.581199354001606143928050809e21,
+       -.6672106568924916298020941484e20,
+       0.2316433580634002297931815435e19,
+       -.3588817569910106050743641413e17,
+       0.2908795263834775409737601689e15,
+       -.1322983480332126453125473247e13,
+       0.3413234182301700539091292655e10,
+       -.4695753530642995859767162166e7,
+       0.2701122710892323414856790990e4,
+};
+static double q1[] = {
+       0.1162398708003212287858529400e22,
+       0.1185770712190320999837113348e20,
+       0.6092061398917521746105196863e17,
+       0.2081661221307607351240184229e15,
+       0.5243710262167649715406728642e12,
+       0.1013863514358673989967045588e10,
+       0.1501793594998585505921097578e7,
+       0.1606931573481487801970916749e4,
+       1.0,
+};
+static double p2[] = {
+       -.4435757816794127857114720794e7,
+       -.9942246505077641195658377899e7,
+       -.6603373248364939109255245434e7,
+       -.1523529351181137383255105722e7,
+       -.1098240554345934672737413139e6,
+       -.1611616644324610116477412898e4,
+       0.0,
+};
+static double q2[] = {
+       -.4435757816794127856828016962e7,
+       -.9934124389934585658967556309e7,
+       -.6585339479723087072826915069e7,
+       -.1511809506634160881644546358e7,
+       -.1072638599110382011903063867e6,
+       -.1455009440190496182453565068e4,
+       1.0,
+};
+static double p3[] = {
+       0.3322091340985722351859704442e5,
+       0.8514516067533570196555001171e5,
+       0.6617883658127083517939992166e5,
+       0.1849426287322386679652009819e5,
+       0.1706375429020768002061283546e4,
+       0.3526513384663603218592175580e2,
+       0.0,
+};
+static double q3[] = {
+       0.7087128194102874357377502472e6,
+       0.1819458042243997298924553839e7,
+       0.1419460669603720892855755253e7,
+       0.4002944358226697511708610813e6,
+       0.3789022974577220264142952256e5,
+       0.8638367769604990967475517183e3,
+       1.0,
+};
+static double p4[] = {
+       -.9963753424306922225996744354e23,
+       0.2655473831434854326894248968e23,
+       -.1212297555414509577913561535e22,
+       0.2193107339917797592111427556e20,
+       -.1965887462722140658820322248e18,
+       0.9569930239921683481121552788e15,
+       -.2580681702194450950541426399e13,
+       0.3639488548124002058278999428e10,
+       -.2108847540133123652824139923e7,
+       0.0,
+};
+static double q4[] = {
+       0.5082067366941243245314424152e24,
+       0.5435310377188854170800653097e22,
+       0.2954987935897148674290758119e20,
+       0.1082258259408819552553850180e18,
+       0.2976632125647276729292742282e15,
+       0.6465340881265275571961681500e12,
+       0.1128686837169442121732366891e10,
+       0.1563282754899580604737366452e7,
+       0.1612361029677000859332072312e4,
+       1.0,
+};
+
+double
+j1(arg) double arg;{
+       double xsq, n, d, x;
+       double sin(), cos(), sqrt();
+       int i;
+
+       x = arg;
+       if(x < 0.) x = -x;
+       if(x > 8.){
+               asympt(x);
+               n = x - 3.*pio4;
+               n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
+               if(arg <0.) n = -n;
+               return(n);
+       }
+       xsq = x*x;
+       for(n=0,d=0,i=8;i>=0;i--){
+               n = n*xsq + p1[i];
+               d = d*xsq + q1[i];
+       }
+       return(arg*n/d);
+}
+
+double
+y1(arg) double arg;{
+       double xsq, n, d, x;
+       double sin(), cos(), sqrt(), log(), j1();
+       int i;
+
+       errno = 0;
+       x = arg;
+       if(x <= 0.){
+               errno = EDOM;
+               return(-HUGE);
+       }
+       if(x > 8.){
+               asympt(x);
+               n = x - 3*pio4;
+               return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
+       }
+       xsq = x*x;
+       for(n=0,d=0,i=9;i>=0;i--){
+               n = n*xsq + p4[i];
+               d = d*xsq + q4[i];
+       }
+       return(x*n/d + tpi*(j1(x)*log(x)-1./x));
+}
+
+static
+asympt(arg) double arg;{
+       double zsq, n, d;
+       int i;
+       zsq = 64./(arg*arg);
+       for(n=0,d=0,i=6;i>=0;i--){
+               n = n*zsq + p2[i];
+               d = d*zsq + q2[i];
+       }
+       pzero = n/d;
+       for(n=0,d=0,i=6;i>=0;i--){
+               n = n*zsq + p3[i];
+               d = d*zsq + q3[i];
+       }
+       qzero = (8./arg)*(n/d);
+}

diff --git a/usr/src/lib/libm/jn.c b/usr/src/lib/libm/jn.c
new file mode 100644 (file)
index 0000000..f9713a0
--- /dev/null
+++ b/usr/src/lib/libm/jn.c
@@ -0,0 +1,107 @@
+/*
+       floating point Bessel's function of
+       the first and second kinds and of
+       integer order.
+
+       int n;
+       double x;
+       jn(n,x);
+
+       returns the value of Jn(x) for all
+       integer values of n and all real values
+       of x.
+
+       There are no error returns.
+       Calls j0, j1.
+
+       For n=0, j0(x) is called,
+       for n=1, j1(x) is called,
+       for n<x, forward recursion us used starting
+       from values of j0(x) and j1(x).
+       for n>x, a continued fraction approximation to
+       j(n,x)/j(n-1,x) is evaluated and then backward
+       recursion is used starting from a supposed value
+       for j(n,x). The resulting value of j(0,x) is
+       compared with the actual value to correct the
+       supposed value of j(n,x).
+
+       yn(n,x) is similar in all respects, except
+       that forward recursion is used for all
+       values of n>1.
+*/
+
+#include <math.h>
+#include <errno.h>
+
+int    errno;
+
+double
+jn(n,x) int n; double x;{
+       int i;
+       double a, b, temp;
+       double xsq, t;
+       double j0(), j1();
+
+       if(n<0){
+               n = -n;
+               x = -x;
+       }
+       if(n==0) return(j0(x));
+       if(n==1) return(j1(x));
+       if(x == 0.) return(0.);
+       if(n>x) goto recurs;
+
+       a = j0(x);
+       b = j1(x);
+       for(i=1;i<n;i++){
+               temp = b;
+               b = (2.*i/x)*b - a;
+               a = temp;
+       }
+       return(b);
+
+recurs:
+       xsq = x*x;
+       for(t=0,i=n+16;i>n;i--){
+               t = xsq/(2.*i - t);
+       }
+       t = x/(2.*n-t);
+
+       a = t;
+       b = 1;
+       for(i=n-1;i>0;i--){
+               temp = b;
+               b = (2.*i/x)*b - a;
+               a = temp;
+       }
+       return(t*j0(x)/b);
+}
+
+double
+yn(n,x) int n; double x;{
+       int i;
+       int sign;
+       double a, b, temp;
+       double y0(), y1();
+
+       if (x <= 0) {
+               errno = EDOM;
+               return(-HUGE);
+       }
+       sign = 1;
+       if(n<0){
+               n = -n;
+               if(n%2 == 1) sign = -1;
+       }
+       if(n==0) return(y0(x));
+       if(n==1) return(sign*y1(x));
+
+       a = y0(x);
+       b = y1(x);
+       for(i=1;i<n;i++){
+               temp = b;
+               b = (2.*i/x)*b - a;
+               a = temp;
+       }
+       return(sign*b);
+}

diff --git a/usr/src/lib/libm/log.c b/usr/src/lib/libm/log.c
new file mode 100644 (file)
index 0000000..3c9d65d
--- /dev/null
+++ b/usr/src/lib/libm/log.c
@@ -0,0 +1,62 @@
+/*
+       log returns the natural logarithm of its floating
+       point argument.
+
+       The coefficients are #2705 from Hart & Cheney. (19.38D)
+
+       It calls frexp.
+*/
+
+#include <errno.h>
+#include <math.h>
+
+int    errno;
+double frexp();
+static double  log2    = 0.693147180559945309e0;
+static double  ln10    = 2.302585092994045684;
+static double  sqrto2  = 0.707106781186547524e0;
+static double  p0      = -.240139179559210510e2;
+static double  p1      = 0.309572928215376501e2;
+static double  p2      = -.963769093368686593e1;
+static double  p3      = 0.421087371217979714e0;
+static double  q0      = -.120069589779605255e2;
+static double  q1      = 0.194809660700889731e2;
+static double  q2      = -.891110902798312337e1;
+
+double
+log(arg)
+double arg;
+{
+       double x,z, zsq, temp;
+       int exp;
+
+       if(arg <= 0.) {
+               errno = EDOM;
+               return(-HUGE);
+       }
+       x = frexp(arg,&exp);
+       while(x<0.5) {
+               x = x*2;
+               exp = exp-1;
+       }
+       if(x<sqrto2) {
+               x = 2*x;
+               exp = exp-1;
+       }
+
+       z = (x-1)/(x+1);
+       zsq = z*z;
+
+       temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0;
+       temp = temp/(((1.0*zsq + q2)*zsq + q1)*zsq + q0);
+       temp = temp*z + exp*log2;
+       return(temp);
+}
+
+double
+log10(arg)
+double arg;
+{
+
+       return(log(arg)/ln10);
+}

diff --git a/usr/src/lib/libm/pow.c b/usr/src/lib/libm/pow.c
new file mode 100644 (file)
index 0000000..b2708db
--- /dev/null
+++ b/usr/src/lib/libm/pow.c
@@ -0,0 +1,36 @@
+/*
+       computes a^b.
+       uses log and exp
+*/
+
+#include       <errno.h>
+int errno;
+double log(), exp();
+
+double
+pow(arg1,arg2)
+double arg1, arg2;
+{
+       double temp;
+       long l;
+
+       if(arg1 <= 0.) {
+               if(arg1 == 0.) {
+                       if(arg2 <= 0.)
+                               goto domain;
+                       return(0.);
+               }
+               l = arg2;
+               if(l != arg2)
+                       goto domain;
+               temp = exp(arg2 * log(-arg1));
+               if(l & 1)
+                       temp = -temp;
+               return(temp);
+       }
+       return(exp(arg2 * log(arg1)));
+
+domain:
+       errno = EDOM;
+       return(0.);
+}

diff --git a/usr/src/lib/libm/sin.c b/usr/src/lib/libm/sin.c
new file mode 100644 (file)
index 0000000..3f76374
--- /dev/null
+++ b/usr/src/lib/libm/sin.c
@@ -0,0 +1,74 @@
+/*
+       C program for floating point sin/cos.
+       Calls modf.
+       There are no error exits.
+       Coefficients are #3370 from Hart & Cheney (18.80D).
+*/
+
+static double twoopi   = 0.63661977236758134308;
+static double p0       =  .1357884097877375669092680e8;
+static double p1       = -.4942908100902844161158627e7;
+static double p2       =  .4401030535375266501944918e6;
+static double p3       = -.1384727249982452873054457e5;
+static double p4       =  .1459688406665768722226959e3;
+static double q0       =  .8644558652922534429915149e7;
+static double q1       =  .4081792252343299749395779e6;
+static double q2       =  .9463096101538208180571257e4;
+static double q3       =  .1326534908786136358911494e3;
+
+double
+cos(arg)
+double arg;
+{
+       double sinus();
+       if(arg<0)
+               arg = -arg;
+       return(sinus(arg, 1));
+}
+
+double
+sin(arg)
+double arg;
+{
+       double sinus();
+       return(sinus(arg, 0));
+}
+
+static double
+sinus(arg, quad)
+double arg;
+int quad;
+{
+       double modf();
+       double e, f;
+       double ysq;
+       double x,y;
+       int k;
+       double temp1, temp2;
+
+       x = arg;
+       if(x<0) {
+               x = -x;
+               quad = quad + 2;
+       }
+       x = x*twoopi;   /*underflow?*/
+       if(x>32764){
+               y = modf(x,&e);
+               e = e + quad;
+               modf(0.25*e,&f);
+               quad = e - 4*f;
+       }else{
+               k = x;
+               y = x - k;
+               quad = (quad + k) & 03;
+       }
+       if (quad & 01)
+               y = 1-y;
+       if(quad > 1)
+               y = -y;
+
+       ysq = y*y;
+       temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y;
+       temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0);
+       return(temp1/temp2);
+}

diff --git a/usr/src/lib/libm/sinh.c b/usr/src/lib/libm/sinh.c
new file mode 100644 (file)
index 0000000..e0e8862
--- /dev/null
+++ b/usr/src/lib/libm/sinh.c
@@ -0,0 +1,67 @@
+/*
+       sinh(arg) returns the hyperbolic sine of its floating-
+       point argument.
+
+       The exponential function is called for arguments
+       greater in magnitude than 0.5.
+
+       A series is used for arguments smaller in magnitude than 0.5.
+       The coefficients are #2029 from Hart & Cheney. (20.36D)
+
+       cosh(arg) is computed from the exponential function for
+       all arguments.
+*/
+
+double exp();
+
+static double p0  = -0.6307673640497716991184787251e+6;
+static double p1  = -0.8991272022039509355398013511e+5;
+static double p2  = -0.2894211355989563807284660366e+4;
+static double p3  = -0.2630563213397497062819489e+2;
+static double q0  = -0.6307673640497716991212077277e+6;
+static double q1   = 0.1521517378790019070696485176e+5;
+static double q2  = -0.173678953558233699533450911e+3;
+
+double
+sinh(arg)
+double arg;
+{
+       double temp, argsq;
+       register sign;
+
+       sign = 1;
+       if(arg < 0) {
+               arg = - arg;
+               sign = -1;
+       }
+
+       if(arg > 21.) {
+               temp = exp(arg)/2;
+               if (sign>0)
+                       return(temp);
+               else
+                       return(-temp);
+       }
+
+       if(arg > 0.5) {
+               return(sign*(exp(arg) - exp(-arg))/2);
+       }
+
+       argsq = arg*arg;
+       temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg;
+       temp /= (((argsq+q2)*argsq+q1)*argsq+q0);
+       return(sign*temp);
+}
+
+double
+cosh(arg)
+double arg;
+{
+       if(arg < 0)
+               arg = - arg;
+       if(arg > 21.) {
+               return(exp(arg)/2);
+       }
+
+       return((exp(arg) + exp(-arg))/2);
+}

diff --git a/usr/src/lib/libm/sqrt.c b/usr/src/lib/libm/sqrt.c
new file mode 100644 (file)
index 0000000..0282333
--- /dev/null
+++ b/usr/src/lib/libm/sqrt.c
@@ -0,0 +1,56 @@
+/*
+       sqrt returns the square root of its floating
+       point argument. Newton's method.
+
+       calls frexp
+*/
+
+#include <errno.h>
+
+int errno;
+double frexp();
+
+double
+sqrt(arg)
+double arg;
+{
+       double x, temp;
+       int exp;
+       int i;
+
+       if(arg <= 0.) {
+               if(arg < 0.)
+                       errno = EDOM;
+               return(0.);
+       }
+       x = frexp(arg,&exp);
+       while(x < 0.5) {
+               x *= 2;
+               exp--;
+       }
+       /*
+        * NOTE
+        * this wont work on 1's comp
+        */
+       if(exp & 1) {
+               x *= 2;
+               exp--;
+       }
+       temp = 0.5*(1.0+x);
+
+       while(exp > 60) {
+               temp *= (1L<<30);
+               exp -= 60;
+       }
+       while(exp < -60) {
+               temp /= (1L<<30);
+               exp += 60;
+       }
+       if(exp >= 0)
+               temp *= 1L << (exp/2);
+       else
+               temp /= 1L << (-exp/2);
+       for(i=0; i<=4; i++)
+               temp = 0.5*(temp + arg/temp);
+       return(temp);
+}

diff --git a/usr/src/lib/libm/tan.c b/usr/src/lib/libm/tan.c
new file mode 100644 (file)
index 0000000..7df1719
--- /dev/null
+++ b/usr/src/lib/libm/tan.c
@@ -0,0 +1,73 @@
+/*
+       floating point tangent
+
+       A series is used after range reduction.
+       Coefficients are #4285 from Hart & Cheney. (19.74D)
+*/
+
+#include <errno.h>
+#include <math.h>
+
+int    errno;
+static double invpi      = 1.27323954473516268;
+static double p0        = -0.1306820264754825668269611177e+5;
+static double p1         = 0.1055970901714953193602353981e+4;
+static double p2        = -0.1550685653483266376941705728e+2;
+static double p3         = 0.3422554387241003435328470489e-1;
+static double p4         = 0.3386638642677172096076369e-4;
+static double q0        = -0.1663895238947119001851464661e+5;
+static double q1         = 0.4765751362916483698926655581e+4;
+static double q2        = -0.1555033164031709966900124574e+3;
+
+double
+tan(arg)
+double arg;
+{
+       double modf();
+       double sign, temp, e, x, xsq;
+       int flag, i;
+
+       flag = 0;
+       sign = 1.;
+       if(arg < 0.){
+               arg = -arg;
+               sign = -1.;
+       }
+       arg = arg*invpi;   /*overflow?*/
+       x = modf(arg,&e);
+       i = e;
+       switch(i%4) {
+       case 1:
+               x = 1. - x;
+               flag = 1;
+               break;
+
+       case 2:
+               sign = - sign;
+               flag = 1;
+               break;
+
+       case 3:
+               x = 1. - x;
+               sign = - sign;
+               break;
+
+       case 0:
+               break;
+       }
+
+       xsq = x*x;
+       temp = ((((p4*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0)*x;
+       temp = temp/(((1.0*xsq+q2)*xsq+q1)*xsq+q0);
+
+       if(flag == 1) {
+               if(temp == 0.) {
+                       errno = ERANGE;
+                       if (sign>0)
+                               return(HUGE);
+                       return(-HUGE);
+               }
+               temp = 1./temp;
+       }
+       return(sign*temp);
+}

diff --git a/usr/src/lib/libm/tanh.c b/usr/src/lib/libm/tanh.c
new file mode 100644 (file)
index 0000000..3f83a3f
--- /dev/null
+++ b/usr/src/lib/libm/tanh.c
@@ -0,0 +1,27 @@
+/*
+       tanh(arg) computes the hyperbolic tangent of its floating
+       point argument.
+
+       sinh and cosh are called except for large arguments, which
+       would cause overflow improperly.
+*/
+
+double sinh(), cosh();
+
+double
+tanh(arg)
+double arg;
+{
+       double sign;
+
+       sign = 1.;
+       if(arg < 0.){
+               arg = -arg;
+               sign = -1.;
+       }
+
+       if(arg > 21.)
+               return(sign);
+
+       return(sign*sinh(arg)/cosh(arg));
+}

diff --git a/usr/src/lib/libmp/gcd.c b/usr/src/lib/libmp/gcd.c
new file mode 100644 (file)
index 0000000..cedad0a
--- /dev/null
+++ b/usr/src/lib/libmp/gcd.c
@@ -0,0 +1,46 @@
+#include <mp.h>
+gcd(a,b,c) MINT *a,*b,*c;
+{      MINT x,y,z,w;
+       x.len=y.len=z.len=w.len=0;
+       move(a,&x);
+       move(b,&y);
+       while(y.len!=0)
+       {       mdiv(&x,&y,&w,&z);
+               move(&y,&x);
+               move(&z,&y);
+       }
+       move(&x,c);
+       xfree(&x);
+       xfree(&y);
+       xfree(&z);
+       xfree(&w);
+       return;
+}
+invert(a, b, c) MINT *a, *b, *c;
+{      MINT x, y, z, w, Anew, Aold;
+       int i = 0;
+       x.len = y.len = z.len = w.len = Aold.len = 0;
+       Anew.len = 1;
+       Anew.val = xalloc(1);
+       *Anew.val = 1;
+       move(b, &x);
+       move(a, &y);
+       while(y.len != 0)
+       {       mdiv(&x, &y, &w, &z);
+               move(&Anew, &x);
+               mult(&w, &Anew, &Anew);
+               madd(&Anew, &Aold, &Anew);
+               move(&x, &Aold);
+               move(&y, &x);
+               move(&z, &y);
+               i++;
+       }
+       move(&Aold, c);
+       if( (i&01) == 0) msub(b, c, c);
+       xfree(&x);
+       xfree(&y);
+       xfree(&z);
+       xfree(&w);
+       xfree(&Aold);
+       xfree(&Anew);
+}

diff --git a/usr/src/lib/libmp/madd.c b/usr/src/lib/libmp/madd.c
new file mode 100644 (file)
index 0000000..155d8c0
--- /dev/null
+++ b/usr/src/lib/libmp/madd.c
@@ -0,0 +1,146 @@
+#include <mp.h>
+m_add(a,b,c) struct mint *a,*b,*c;
+{      int carry,i;
+       int x;
+       short *cval;
+       cval=xalloc(a->len+1,"m_add");
+       carry=0;
+       for(i=0;i<b->len;i++)
+       {       x=carry+a->val[i]+b->val[i];
+               if(x&0100000)
+               {       carry=1;
+                       cval[i]=x&077777;
+               }
+               else
+               {       carry=0;
+                       cval[i]=x;
+               }
+       }
+       for(;i<a->len;i++)
+       {       x=carry+a->val[i];
+               if(x&0100000) cval[i]=x&077777;
+               else
+               {       carry=0;
+                       cval[i]=x;
+               }
+
+       }
+       if(carry==1)
+       {       cval[i]=1;
+               c->len=i+1;
+       }
+       else c->len=a->len;
+       c->val=cval;
+       if(c->len==0) shfree(cval);
+       return;
+}
+madd(a,b,c) struct mint *a,*b,*c;
+{      struct mint x,y,z;
+       int sign;
+       x.len=a->len;
+       x.val=a->val;
+       y.len=b->len;
+       y.val=b->val;
+       z.len=0;
+       sign=1;
+       if(x.len>=0)
+               if(y.len>=0)
+                       if(x.len>=y.len) m_add(&x,&y,&z);
+                       else m_add(&y,&x,&z);
+               else
+               {       y.len= -y.len;
+                       msub(&x,&y,&z);
+               }
+       else    if(y.len<=0)
+               {       x.len = -x.len;
+                       y.len= -y.len;
+                       sign= -1;
+                       madd(&x,&y,&z);
+               }
+               else
+               {       x.len= -x.len;
+                       msub(&y,&x,&z);
+               }
+        xfree(c);
+       c->val=z.val;
+       c->len=sign*z.len;
+       return;
+}
+m_sub(a,b,c) struct mint *a,*b,*c;
+{      int x,i;
+       int borrow;
+       short one;
+       struct mint mone;
+       one=1; mone.len= 1; mone.val= &one;
+       c->val=xalloc(a->len,"m_sub");
+       borrow=0;
+       for(i=0;i<b->len;i++)
+       {       x=borrow+a->val[i]-b->val[i];
+               if(x&0100000)
+               {       borrow= -1;
+                       c->val[i]=x&077777;
+               }
+               else
+               {       borrow=0;
+                       c->val[i]=x;
+               }
+       }
+       for(;i<a->len;i++)
+       {       x=borrow+a->val[i];
+               if(x&0100000) c->val[i]=x&077777;
+               else
+               {       borrow=0;
+                       c->val[i]=x;
+               }
+       }
+       if(borrow<0)
+       {       for(i=0;i<a->len;i++) c->val[i] ^= 077777;
+               c->len=a->len;
+               madd(c,&mone,c);
+       }
+       for(i=a->len-1;i>=0;--i) if(c->val[i]>0)
+                               {       if(borrow==0) c->len=i+1;
+                                       else c->len= -i-1;
+                                       return;
+                               }
+       shfree(c->val);
+       return;
+}
+msub(a,b,c) struct mint *a,*b,*c;
+{      struct mint x,y,z;
+       int sign;
+       x.len=a->len;
+       y.len=b->len;
+       x.val=a->val;
+       y.val=b->val;
+       z.len=0;
+       sign=1;
+       if(x.len>=0)
+               if(y.len>=0)
+                       if(x.len>=y.len) m_sub(&x,&y,&z);
+                       else
+                       {       sign= -1;
+                               msub(&y,&x,&z);
+                       }
+               else
+               {       y.len= -y.len;
+                       madd(&x,&y,&z);
+               }
+       else    if(y.len<=0)
+               {       sign= -1;
+                       x.len= -x.len;
+                       y.len= -y.len;
+                       msub(&y,&x,&z);
+               }
+               else
+               {       x.len= -x.len;
+                       madd(&x,&y,&z);
+                       sign= -1;
+               }
+       if(a==c && x.len!=0) xfree(a);
+       else if(b==c && y.len!=0) xfree(b);
+       else xfree(c);
+       c->val=z.val;
+       c->len=sign*z.len;
+       return;
+}

diff --git a/usr/src/lib/libmp/mdiv.c b/usr/src/lib/libmp/mdiv.c
new file mode 100644 (file)
index 0000000..1f8a0fc
--- /dev/null
+++ b/usr/src/lib/libmp/mdiv.c
@@ -0,0 +1,106 @@
+#include <mp.h>
+mdiv(a,b,q,r) MINT *a,*b,*q,*r;
+{      MINT x,y;
+       int sign;
+       sign=1;
+       x.val=a->val;
+       y.val=b->val;
+       x.len=a->len;
+       if(x.len<0) {sign= -1; x.len= -x.len;}
+       y.len=b->len;
+       if(y.len<0) {sign= -sign; y.len= -y.len;}
+       xfree(q);
+       xfree(r);
+       m_div(&x,&y,q,r);
+       if(sign==-1)
+       {       q->len= -q->len;
+               r->len = - r->len;
+       }
+       return;
+}
+m_dsb(q,n,a,b) short *a,*b;
+{      long int x,qx;
+       int borrow,j,u;
+       qx=q;
+       borrow=0;
+       for(j=0;j<n;j++)
+       {       x=borrow-a[j]*qx+b[j];
+               b[j]=x&077777;
+               borrow=x>>15;
+       }
+       x=borrow+b[j];
+       b[j]=x&077777;
+       if(x>>15 ==0) { return(0);}
+       borrow=0;
+       for(j=0;j<n;j++)
+       {       u=a[j]+b[j]+borrow;
+               if(u<0) borrow=1;
+               else borrow=0;
+               b[j]=u&077777;
+       }
+       { return(1);}
+}
+m_trq(v1,v2,u1,u2,u3)
+{      long int d;
+       long int x1;
+       if(u1==v1) d=077777;
+       else d=(u1*0100000L+u2)/v1;
+       while(1)
+       {       x1=u1*0100000L+u2-v1*d;
+               x1=x1*0100000L+u3-v2*d;
+               if(x1<0) d=d-1;
+               else {return(d);}
+       }
+}
+m_div(a,b,q,r) MINT *a,*b,*q,*r;
+{      MINT u,v,x,w;
+       short d,*qval;
+       int qq,n,v1,v2,j;
+       u.len=v.len=x.len=w.len=0;
+       if(b->len==0) { fatal("mdiv divide by zero"); return;}
+       if(b->len==1)
+       {       r->val=xalloc(1,"m_div1");
+               sdiv(a,b->val[0],q,r->val);
+               if(r->val[0]==0)
+               {       shfree(r->val);
+                       r->len=0;
+               }
+               else r->len=1;
+               return;
+       }
+       if(a->len < b->len)
+       {       q->len=0;
+               r->len=a->len;
+               r->val=xalloc(r->len,"m_div2");
+               for(qq=0;qq<r->len;qq++) r->val[qq]=a->val[qq];
+               return;
+       }
+       x.len=1;
+       x.val = &d;
+       n=b->len;
+       d=0100000L/(b->val[n-1]+1L);
+       mult(a,&x,&u); /*subtle: relies on fact that mult allocates extra space */
+       mult(b,&x,&v);
+       v1=v.val[n-1];
+       v2=v.val[n-2];
+       qval=xalloc(a->len-n+1,"m_div3");
+       for(j=a->len-n;j>=0;j--)
+       {       qq=m_trq(v1,v2,u.val[j+n],u.val[j+n-1],u.val[j+n-2]);
+               if(m_dsb(qq,n,v.val,&(u.val[j]))) qq -= 1;
+               qval[j]=qq;
+       }
+       x.len=n;
+       x.val=u.val;
+       mcan(&x);
+       sdiv(&x,d,&w,(short *)&qq);
+       r->len=w.len;
+       r->val=w.val;
+       q->val=qval;
+       qq=a->len-n+1;
+       if(qq>0 && qval[qq-1]==0) qq -= 1;
+       q->len=qq;
+       if(qq==0) shfree(qval);
+       if(x.len!=0) xfree(&u);
+       xfree(&v);
+       return;
+}

diff --git a/usr/src/lib/libmp/mout.c b/usr/src/lib/libmp/mout.c
new file mode 100644 (file)
index 0000000..afb30a8
--- /dev/null
+++ b/usr/src/lib/libmp/mout.c
@@ -0,0 +1,141 @@
+#include <stdio.h>
+#include <mp.h>
+m_in(a,b,f) MINT *a; FILE *f;
+{      MINT x,y,ten;
+       int sign,c;
+       short qten,qy;
+       xfree(a);
+       sign=1;
+       ten.len=1;
+       ten.val= &qten;
+       qten=b;
+       x.len=0;
+       y.len=1;
+       y.val= &qy;
+       while((c=getc(f))!=EOF)
+       switch(c)
+       {
+       case '\\':      getc(f);
+               continue;
+       case '\t':
+       case '\n': a->len *= sign;
+               xfree(&x);
+               return(0);
+       case ' ':
+               continue;
+       case '-': sign = -sign;
+               continue;
+       default: if(c>='0' && c<= '9')
+               {       qy=c-'0';
+                       mult(&x,&ten,a);
+                       madd(a,&y,a);
+                       move(a,&x);
+                       continue;
+               }
+               else
+               {       VOID ungetc(c,stdin);
+                       a->len *= sign;
+                       return(0);
+               }
+       }
+       return(EOF);
+}
+m_out(a,b,f) MINT *a; FILE *f;
+{      int sign,xlen,i;
+       short r;
+       MINT x;
+       char *obuf;
+       register char *bp;
+       sign=1;
+       xlen=a->len;
+       if(xlen<0)
+       {       xlen= -xlen;
+               sign= -1;
+       }
+       if(xlen==0)
+       {       fprintf(f,"0\n");
+               return;
+       }
+       x.len=xlen;
+       x.val=xalloc(xlen,"m_out");
+       for(i=0;i<xlen;i++) x.val[i]=a->val[i];
+       obuf=(char *)malloc(7*xlen);
+       bp=obuf+7*xlen-1;
+       *bp--=0;
+       while(x.len>0)
+       {       for(i=0;i<10&&x.len>0;i++)
+               {       sdiv(&x,b,&x,&r);
+                       *bp--=r+'0';
+               }
+               if(x.len>0) *bp--=' ';
+       }
+       if(sign==-1) *bp--='-';
+       fprintf(f,"%s\n",bp+1);
+       free(obuf);
+       FREE(x)
+       return;
+}
+sdiv(a,n,q,r) MINT *a,*q; short *r;
+{      MINT x,y;
+       int sign;
+       sign=1;
+       x.len=a->len;
+       x.val=a->val;
+       if(n<0)
+       {       sign= -sign;
+               n= -n;
+       }
+       if(x.len<0)
+       {       sign = -sign;
+               x.len= -x.len;
+       }
+       s_div(&x,n,&y,r);
+       xfree(q);
+       q->val=y.val;
+       q->len=sign*y.len;
+       *r = *r*sign;
+       return;
+}
+s_div(a,n,q,r) MINT *a,*q; short *r;
+{      int qlen,i;
+       long int x;
+       short *qval;
+       x=0;
+       qlen=a->len;
+       qval=xalloc(qlen,"s_div");
+       for(i=qlen-1;i>=0;i--)
+       {
+               x=x*0100000L+a->val[i];
+               qval[i]=x/n;
+               x=x%n;
+       }
+       *r=x;
+       if(qval[qlen-1]==0) qlen--;
+       q->len=qlen;
+       q->val=qval;
+       if(qlen==0) shfree(qval);
+       return;
+}
+min(a) MINT *a;
+{
+       return(m_in(a,10,stdin));
+}
+omin(a) MINT *a;
+{
+       return(m_in(a,8,stdin));
+}
+mout(a) MINT *a;
+{
+       m_out(a,10,stdout);
+}
+omout(a) MINT *a;
+{
+       m_out(a,8,stdout);
+}
+fmout(a,f) MINT *a; FILE *f;
+{      m_out(a,10,f);
+}
+fmin(a,f) MINT *a; FILE *f;
+{
+       return(m_in(a,10,f));
+}

diff --git a/usr/src/lib/libmp/msqrt.c b/usr/src/lib/libmp/msqrt.c
new file mode 100644 (file)
index 0000000..a954593
--- /dev/null
+++ b/usr/src/lib/libmp/msqrt.c
@@ -0,0 +1,37 @@
+#include <mp.h>
+msqrt(a,b,r) MINT *a,*b,*r;
+{      MINT x,junk,y;
+       int j;
+       x.len=junk.len=y.len=0;
+       if(a->len<0) fatal("msqrt: neg arg");
+       if(a->len==0)
+       {       b->len=0;
+               r->len=0;
+               return(0);
+       }
+       if(a->len%2==1) x.len=(1+a->len)/2;
+       else x.len=1+a->len/2;
+       x.val=xalloc(x.len,"msqrt");
+       for(j=0;j<x.len;x.val[j++]=0);
+       if(a->len%2==1) x.val[x.len-1]=0400;
+       else x.val[x.len-1]=1;
+       xfree(b);
+       xfree(r);
+loop:
+       mdiv(a,&x,&y,&junk);
+       xfree(&junk);
+       madd(&x,&y,&y);
+       sdiv(&y,2,&y,(short *)&j);
+       if(mcmp(&x,&y)>0)
+       {       xfree(&x);
+               move(&y,&x);
+               xfree(&y);
+               goto loop;
+       }
+       xfree(&y);
+       move(&x,b);
+       mult(&x,&x,&x);
+       msub(a,&x,r);
+       xfree(&x);
+       return(r->len);
+}

diff --git a/usr/src/lib/libmp/mult.c b/usr/src/lib/libmp/mult.c
new file mode 100644 (file)
index 0000000..d4d8de9
--- /dev/null
+++ b/usr/src/lib/libmp/mult.c
@@ -0,0 +1,87 @@
+#include <mp.h>
+mult(a,b,c) struct mint *a,*b,*c;
+{      struct mint x,y,z;
+       int sign;
+       sign = 1;
+       x.val=a->val;
+       y.val=b->val;
+       z.len=0;
+       if(a->len<0)
+       {       x.len= -a->len;
+               sign= -sign;
+       }
+       else    x.len=a->len;
+       if(b->len<0)
+       {       y.len= -b->len;
+               sign= -sign;
+       }
+       else    y.len=b->len;
+       if(x.len<y.len) m_mult(&y,&x,&z);
+       else m_mult(&x,&y,&z);
+       xfree(c);
+       if(sign<0) c->len= -z.len;
+       else c->len=z.len;
+       if(c->len==0) shfree(z.val);
+       else c->val=z.val;
+       return;
+}
+#define S2 x=a->val[j];
+#define S3 x=x*b->val[i-j];
+#define S4 tradd(&carry,&sum,x);
+#define S5 c->val[i]=sum.yy.low&077777;
+#define S6 sum.xx=sum.xx>>15;
+#define S7 sum.yy.high=carry;
+m_mult(a,b,c) struct mint *a,*b,*c;
+{      long x;
+       union {long xx; struct half yy;} sum;
+       int carry;
+       int i,j;
+       c->val=xalloc(a->len+b->len,"m_mult");
+       sum.xx=0;
+       for(i=0;i<b->len;i++)
+       {       carry=0;
+               for(j=0;j<i+1;j++)
+               {       S2
+                       S3
+                       S4
+               }
+               S5
+               S6
+               S7
+       }
+       for(;i<a->len;i++)
+       {       carry=0;
+               for(j=i-b->len+1;j<i+1;j++)
+               {       S2
+                       S3
+                       S4
+               }
+               S5
+               S6
+               S7
+       }
+       for(;i<a->len+b->len;i++)
+       {       carry=0;
+               for(j=i-b->len+1;j<a->len;j++)
+               {       S2
+                       S3
+                       S4
+               }
+               S5
+               S6
+               S7
+       }
+       if(c->val[i-1]!=0)
+               c->len=a->len+b->len;
+       else    c->len=a->len+b->len-1;
+       return;
+}
+tradd(a,b,c) long c; int *a; union g {long xx; struct half yy;} *b;
+{
+       b->xx= b->xx+c;
+       if(b->yy.high&0100000)
+       {       b->yy.high= b->yy.high&077777;
+               *a += 1;
+       }
+       return;
+}

diff --git a/usr/src/lib/libmp/pow.c b/usr/src/lib/libmp/pow.c
new file mode 100644 (file)
index 0000000..4a745e4
--- /dev/null
+++ b/usr/src/lib/libmp/pow.c
@@ -0,0 +1,40 @@
+#include <mp.h>
+pow(a,b,c,d) MINT *a,*b,*c,*d;
+{      int i,j,n;
+       MINT x,y;
+       x.len=y.len=0;
+       xfree(d);
+       d->len=1;
+       d->val=xalloc(1,"pow");
+       *d->val=1;
+       for(j=0;j<b->len;j++)
+       {       n=b->val[b->len-j-1];
+               for(i=0;i<15;i++)
+               {       mult(d,d,&x);
+                       mdiv(&x,c,&y,d);
+                       if((n=n<<1)&0100000)
+                       {       mult(a,d,&x);
+                               mdiv(&x,c,&y,d);
+                       }
+               }
+       }
+       xfree(&x);
+       xfree(&y);
+       return;
+}
+rpow(a,n,b) MINT *a,*b;
+{      MINT x,y;
+       int i;
+       x.len=1;
+       x.val=xalloc(1,"rpow");
+       *x.val=n;
+       y.len=n*a->len+4;
+       y.val=xalloc(y.len,"rpow2");
+       for(i=0;i<y.len;i++) y.val[i]=0;
+       y.val[y.len-1]=010000;
+       xfree(b);
+       pow(a,&x,&y,b);
+       xfree(&x);
+       xfree(&y);
+       return;
+}

diff --git a/usr/src/lib/libmp/util.c b/usr/src/lib/libmp/util.c
new file mode 100644 (file)
index 0000000..6d81d14
--- /dev/null
+++ b/usr/src/lib/libmp/util.c
@@ -0,0 +1,88 @@
+char *malloc();
+#ifdef lint
+int xv_oid;
+#endif
+#include <stdio.h>
+#include <mp.h>
+move(a,b) MINT *a,*b;
+{      int i,j;
+       xfree(b);
+       b->len=a->len;
+       if((i=a->len)<0) i = -i;
+       if(i==0) return;
+       b->val=xalloc(i,"move");
+       for(j=0;j<i;j++)
+               b->val[j]=a->val[j];
+       return;
+}
+dummy(){}
+short *xalloc(nint,s) char *s;
+{      short *i;
+       i=(short *)malloc(2*(unsigned)nint+4);
+#ifdef DBG
+       if(dbg) fprintf(stderr, "%s: %o\n",s,i);
+#endif
+       if(i!=NULL) return(i);
+       fatal("mp: no free space");
+       return(0);
+}
+fatal(s) char *s;
+{
+       fprintf(stderr,"%s\n",s);
+       VOID fflush(stdout);
+       sleep(2);
+       abort();
+}
+xfree(c) MINT *c;
+{
+#ifdef DBG
+       if(dbg) fprintf(stderr, "xfree ");
+#endif
+       if(c->len==0) return;
+       shfree(c->val);
+       c->len=0;
+       return;
+}
+mcan(a) MINT *a;
+{      int i,j;
+       if((i=a->len)==0) return;
+       else if(i<0) i= -i;
+       for(j=i;j>0 && a->val[j-1]==0;j--);
+       if(j==i) return;
+       if(j==0)
+       {       xfree(a);
+               return;
+       }
+       if(a->len > 0) a->len=j;
+       else a->len = -j;
+}
+MINT *itom(n)
+{      MINT *a;
+       a=(MINT *)xalloc(2,"itom");
+       if(n>0)
+       {       a->len=1;
+               a->val=xalloc(1,"itom1");
+               *a->val=n;
+               return(a);
+       }
+       else if(n<0)
+       {       a->len = -1;
+               a->val=xalloc(1,"itom2");
+               *a->val= -n;
+               return(a);
+       }
+       else
+       {       a->len=0;
+               return(a);
+       }
+}
+mcmp(a,b) MINT *a,*b;
+{      MINT c;
+       int res;
+       if(a->len!=b->len) return(a->len-b->len);
+       c.len=0;
+       msub(a,b,&c);
+       res=c.len;
+       xfree(&c);
+       return(res);
+}

diff --git a/usr/src/lib/libnm/pow.c b/usr/src/lib/libnm/pow.c
new file mode 100644 (file)
index 0000000..b2708db
--- /dev/null
+++ b/usr/src/lib/libnm/pow.c
@@ -0,0 +1,36 @@
+/*
+       computes a^b.
+       uses log and exp
+*/
+
+#include       <errno.h>
+int errno;
+double log(), exp();
+
+double
+pow(arg1,arg2)
+double arg1, arg2;
+{
+       double temp;
+       long l;
+
+       if(arg1 <= 0.) {
+               if(arg1 == 0.) {
+                       if(arg2 <= 0.)
+                               goto domain;
+                       return(0.);
+               }
+               l = arg2;
+               if(l != arg2)
+                       goto domain;
+               temp = exp(arg2 * log(-arg1));
+               if(l & 1)
+                       temp = -temp;
+               return(temp);
+       }
+       return(exp(arg2 * log(arg1)));
+
+domain:
+       errno = EDOM;
+       return(0.);
+}