# `#` means: (nothing)

/*
for x from -10 to 10 by 0.5 do
   for y from -10 to 10 by 0.5 do
      lprint(`cosc`, `;`, x,  `+`, y, `*I;`, traperror(evalf(f(x+I*y))))
      #lprint(`cosc`, `;`, x,  `;`, traperror(evalf(f(x))))
   od
od;

Digits := 30:
fd := open(`c:/temp/testdata.txt`, WRITE);
for x from -10 to 10 by 0.5 do
#   for y from -10 to 10 by 0.5 do
      z := x + I*y;
#fprintf(fd, `tanc;%a;%a\n`, z,traperror(evalf(tan(z)/z)))
fprintf(fd, `tanc;%a;%a\n`, x,traperror(evalf(tan(x)/x)))
#   od
od:
close(fd);
*/

mathfns := {
   `convertbase`,
   `gcd`,
   `isprime`,
   `prevprime`,
   `nextprime`,
   `lcm`,
   `max`,
   `min`,
   `Phi`,
   `random`,
   `randomseed`,
   `toradians`,
   `pochhammer` # -> test it independently,put test cases from poch2.tst into math2.tst
}


calcfns := {
   `Si`,
   `Ci`,
   `Shi`,
   `Chi`,
   `Ssi`,
   `dawson`,
   `Ei`,
   `Psi`,
   `dilog`,
   `fresnelc`,
   `fresnels`,
   `xpdiff`,
   `diff`
}


statfns := {
   `ad`,
   `amean`,
   `colnorm`,
   `countentries`,
   `cumsum`,
   `deltalist`,
   `gmean`,
   `hmean`,
   `ios`,
   `iqr`,
   `issorted`,
   `mean`,
   `median`,
   `minmax`,
   `mode`,
   `moment`,
   `numbperm`,
   `obcount`,
   `percentile`,
   `prange`,
   `qmean`,
   `quartiles`,
   `rownorm`,
   `scale`,
   `sd`,
   `skewness`,
   `sorted`,
   `sumdata`,
   `tovals`,
   `var`
}


mathops := {
   `abs`,
   `arccos`,
   `arcsin`,
   `arctan`,
   `cos`,
   `cosh`,
   `entier`,
   `even`,
   `exp`,
   `finite`,
   `lngamma`,
   `int`,
   `imag`,
   `ln`,
   `real`,
   `sign`,
   `sin`,
   `sinh`,
   `sqrt`,
   `tan`,
   `tanh`,
   # baselib functions
   `arccosh`,
   `approx`,
   `arccoth`,
   `arctan2`,
   `arcsinh`,
   `arctanh`,
   `argument`,
   `binomial`,
   `ceil`,
   `conjugate`,
   `cot`,
   `coth`,
   `csc`,
   `fact`,
   `frac`,
   `frexp`,
   `heaviside`,
   `hypot`,
   `irem`,
   `ldexp`,
   `log`,
   `log10`,
   `modf`,
   `round`,
   `sec`,
   `erf`,
   `erfc`,
   `besselj`,
   `bessely`,
   `arcsec`,
   `arccsc`,
   `arccot`,
   `arccsc`,
   `arccsch`,
   `arccot`,
   `arcsech`,
   `csch`,
   `sech`,
   `beta`,
   `root`,
   `proot`,
   `mdf`,
   `xdf`
}

irem := << x, y -> x symmod y >>;  # deprecated but still used in tests

shift := << x, y -> if y > 0 then x <<< y else x >>> y fi >>; # deprecated but still used in tests

mathall := mathfns union mathops;

isdiff := proc(agena, maple, eps) is
   local app;
   if eps = null then
      app := not(approx(agena, maple))
   else
      app := not(approx(agena, maple, eps))
   fi;
   return app and
      not((agena=undefined) and ((maple=undefined) or (type(maple)=complex))) and
      not(agena=infinity and maple=infinity) and not(agena=-infinity and maple=-infinity)
end;

_time := time();

import calc, stats, linalg, cordic, divs;

# test functions with one arg
print("\nChecking functions with one argument\n");

# cosh   ;   0   ;   1.

for keys j in io.lines('math.tst') do
   s := replace(j, ' ', '');  # delete white spaces
   se := s split ';';         # split into its components
   maple := tonumber(se[3]);  # get result
   if se[1] in mathfns then
      prefix := 'math.'
   elif se[1] in calcfns then
      prefix := 'calc.'
   else
      prefix := ''
   fi;
   #prefix := is se[1] in mathfns then 'math.' else '' si;
   s := prefix & se[1] & '(' & se[2] & ')';  # build Agena call
   loadstring('agena := ' & s)();
   if agena = undefined and size se = 3 then  # try complex domain
      s := prefix & se[1] & '(' & se[2] & '+0*I)';  # build Agena call
      loadstring('agena := ' & s)();
   fi
   if maple = 'singularityencountered' or
      maple = 'divisionbyzero' then maple := undefined fi;
   if maple = 'true' then maple := true fi;
   if maple = 'false' then maple := false fi;
#print(s, agena, maple);
   if type(maple) <> boolean then
      #if not(approx(agena, maple, 1e-2)) and not(agena=undefined and maple=undefined) and se[1] <> 'lngamma' then
      if isdiff(agena, maple) then # and se[1] <> 'lngamma' then
         print('>>> ', s, agena, maple)
      fi
   else
      if agena <> maple then
         print('>>> ', s, agena, maple)
      fi
   fi
od;


lb := << x -> ln(x)/ln(2) >>;
for x from -5 to 5 by .1 do
   if not approx(log2(x), lb(x)) then
      print('error in lb, ' & x)
   fi;
   if not approx(log2(x!x), lb(x!x)) then
      print('error in lb, ' & x & '!' & x)
   fi;
od;

for i from -2 to 2 by 0.1 do
   if 1/i <> recip(i) then
      print('Error in recip with real ', i)
   fi
od;

for x from -10 to 10 by 0.1 do
   if not(cis(x) ~= exp(I*x)) then print('Error in cis, ' & x) fi
od;

for x from -10 to 10 by 0.1 do
   z := x + 2*x*I;
   if not(cis(z) ~= cos(z) + I*sin(z)) then print('Error in cis, ' & x) fi
od;

clear z;


# now test functions with two args
print("\nChecking functions with two arguments\n");

# arctan2   ;   1;    0   ;   1.
#Digits := 18:
#fh := open(`c:/p2.txt`, WRITE):
# for x from -10 to 10 by 0.1 do
#    for y from -10 to 10 by 0.1 do
#       fprintf(fh, `%s;%f;%f;%a;no\n`, `pochhammer`, x,  y, traperror(evalf(pochhammer(x, y))))
#    od
# od:
# close(fh);

for keys j in io.lines('math2.tst') do
   s := replace(j, ' ', '');  # delete white spaces
   se := s split ';';         # split into its components
   maple := tonumber(se[4]);  # get result
   if se[1] in mathfns then
      prefix := 'math.'
   elif se[1] in calcfns then
      prefix := 'calc.'
   elif se[1] in statfns then
      prefix := 'stats.'
   else
      prefix := ''
   fi;
   s := prefix & se[1] & '(' & se[2] & ', ' & se[3] & ')';  # build Agena call
   loadstring('agena := ' & s)();
   if agena = undefined and size se = 4 then
      s := prefix & se[1] & '(' & se[2] & '+0*I, ' & se[3] & '+0*I)';  # build Agena call
      loadstring('agena := ' & s)();
   fi;
   if maple = 'singularityencountered' or
      maple = 'divisionbyzero' then maple := undefined fi;
   #if not(approx(agena, maple)) and not((agena=undefined) and (maple=undefined)) then
   if isdiff(agena, maple) then
      print('>>> ', s, agena, maple)
   fi
od;

#/
# now test functions with two args
print("\nChecking functions with three arguments\n");

# diff   ;   1;    0   ;  0   ;   1.

for keys j in io.lines('math3arg.tst') do
   s := replace(j, ' ', '');  # delete white spaces
   se := s split ';';         # split into its components
   maple := tonumber(se[5]);  # get result
   if se[1] in mathfns then
      prefix := 'math.'
   elif se[1] in calcfns then
      prefix := 'calc.'
   else
      prefix := ''
   fi;
   s := prefix & se[1] & '(' & se[2] & ', ' & se[3] & ', ' & se[4] & ')';  # build Agena call
   loadstring('agena := ' & s)();
   if maple = 'singularityencountered' or
      maple = 'divisionbyzero' then maple := undefined fi;
   #if not(approx(agena, maple)) and not((agena=undefined) and (maple=undefined)) then
   if isdiff(agena, maple) then
      print('>>> ', s, agena, maple)
   fi
od;
/#

print("\nChecking Bessel functions\n");

#besselj;0;-10.0-10.0*I;-2314.9753144452130+411.56285702538052*I;no

for keys j in io.lines('bessel.tst') do
   s := replace(j, ' ', '');  # delete white spaces
   se := s split ';';         # split into its components
   maple := tonumber(se[4]);  # get result
   prefix := if se[1] in mathfns then 'math.' else '' fi;
   s := prefix & se[1] & '(' & se[2] & ', ' & se[3] & ')';  # build Agena call
   loadstring('agena := ' & s)();
   if agena = undefined and size se = 4 then
      s := prefix & se[1] & '(' & se[2] & '+0*I, ' & se[3] & '+0*I)';  # build Agena call
      loadstring('agena := ' & s)();
   fi;
   if maple = 'singularityencountered' or
      maple = 'divisionbyzero' then maple := undefined fi;
   if isdiff(agena, maple) then
      # with very large values, roundoff errors are commonplace, so do some dirty comparisons
      if tostring(agena) <> tostring(maple) and (not(agena=undefined and type(maple) = complex)) then
         print('>>> ', s, agena, maple)
      fi
   fi
od;


for i from -2 to 2 by 0.1 do
   if not approx(2^i, antilog2(i)) then
      print('Error in antilog2 with real ', i, 2^i, antilog2(i))
   fi
od;


for i to 100 do
   if ilog2(i) <> entier(log2(i)) then
      print('Error in ilog2 ', i, ilog2(i), entier(log2(i)), ilog2(i) - entier(log2(i)))
   fi
od;


for i to 100 do
   if ilog10(i) <> entier(log10(i)) then
      print('Error in math.ilog10 ', i)
   fi
od;


for i from -100 to 100 do
   if math.copysign(i, i) <> sign(i) * abs(i) then
      print('Error in math.copysign ', i)
   fi
od;

# complex arithmetic: *, -, *, /

print("\nDoing basic complex math\n");

for keys x in io.lines('cmplx1.tst') do
   s := replace(x, ' ', '');  # delete white spaces
   se := s split ';';         # split into its components
   maple := tonumber(se[4]);  # get result
   case se[1]
      of 'a+b' then
         i := tonumber(se[2]);
         j := tonumber(se[3]);
         a := i + I*j;
         b := (i-1) + ((j-1) * I);
         agena := a + b;
         if isdiff(agena, maple) then
            print('>>> ', s, agena, maple)
         fi
      of 'a-b' then
         i := tonumber(se[2]);
         j := tonumber(se[3]);
         a := i + I*j;
         b := (i-1) + ((j-1) * I);
         agena := a - b;
         if isdiff(agena, maple) then
            print('>>> ', s, agena, maple)
         fi
      of 'a*b' then
         i := tonumber(se[2]);
         j := tonumber(se[3]);
         a := i + I*j;
         b := (i-1) + ((j-1) * I);
         agena := a * b;
         if isdiff(agena, maple) then
            print('>>> ', s, agena, maple)
         fi
      of 'a/b' then
         i := tonumber(se[2]);
         j := tonumber(se[3]);
         a := i + I*j;
         b := (i-1) + ((j-1) * I);
         agena := a / b;
         if isdiff(agena, maple) then
            print('>>> ', s, agena, maple)
         fi
      of 'aI^b' then
         i := tonumber(se[2]);
         j := tonumber(se[3]);
         a := i + I*j;
         b := j;
         agena := a ^ b;
         if isdiff(agena, maple) then
            print('>>> ', s, agena, maple)
         fi
      of 'aI^bI' then
         i := tonumber(se[2]);
         j := tonumber(se[3]);
         a := i + I*j;
         b := i + I*j;
         agena := a ^ b;
         if isdiff(agena, maple) then
            print('>>> ', s, agena, maple)
         fi
      of 'a^bI' then
         i := tonumber(se[2]);
         j := tonumber(se[3]);
         a := i;
         b := i + I*j;
         agena := a ^ b;
         if isdiff(agena, maple) then
            print('>>> ', s, agena, maple)
         fi
      of 'aI**j' then
         i := tonumber(se[2]);
         j := tonumber(se[3]);
         a := i + I*j;
         b := j;
         agena := a ** b;
         if isdiff(agena, maple) then
            print('>>> ', s, agena, maple)
         fi
   esac
od;

for i from -2 to 2 by 0.1 do
   for j from -2 to 2 by 0.1 do
      if 1/(i!j) ~<> recip(i!j) then
         print('Error in recip with complex ', i!j)
      fi
   od
od;


for i from -2 to 2 by 0.1 do
   for j from -2 to 2 by 0.1 do
      if not approx(2^(i!j), antilog2(i!j)) then
         print('Error in antilog2 with complex ', i!j)
      fi
   od
od;


for i from -2 to 2 by 0.1 do
   for j from -2 to 2 by 0.1 do
      if not approx(10^(i!j), antilog10(i!j)) then
         print('Error in antilog10 with complex ', i!j)
      fi
   od
od;


print("\nChecking further complex operators\n");
for keys j in io.lines('cmplx2.tst') do
   s := replace(j, ' ', '');  # delete white spaces
   se := s split ';';         # split into its components
   maple := tonumber(se[4]);  # get result
   if se[1] in mathfns then
      prefix := 'math.'
   elif se[1] in calcfns then
      prefix := 'calc.'
   else
      prefix := ''
   fi;
   s := prefix & se[1] & '(' & se[2] & '+I*' & se[3] & ')';  # build Agena call
   loadstring('agena := ' & s)();
   if maple = 'singularityencountered' or
      maple = 'divisionbyzero' then maple := undefined fi;
   #if not(approx(agena, maple)) and not((agena=undefined) and (maple=undefined)) then
   if not((agena=undefined) and (maple=undefined)) and isdiff(agena, maple) then
      print('>>> ', s, agena, maple)
   fi
od;


to 10 do
   if [polar(-1-1*I)] <> [sqrt(2), -3/4*Pi] then
      print('error in polar, 1')
   fi;
   if [polar(1+1*I)] <> [sqrt(2), Pi/4] then
      print('error in polar, 2')
   fi;
   if [polar(1+0*I)] <> [1, 0] then
      print('error in polar, 3')
   fi;
   if [polar(0+0*I)] <> [0, 0] then
      print('error in polar, 4')
   fi
od;


# test functions with one arg
print("\nChecking more elaborate arithmetic in the complex domain\n");

cmplxtest := proc(file::string, f::procedure) is
   local s, se, argu, maple, agena;
   for keys j in io.lines(file) do
      if j[1] = '#' then skip fi;
      s := replace(j, ' ', '');  # delete white spaces
      se := s split ';';         # split into its components
      argu := tonumber(se[1]);   # get argument
      maple := tonumber(se[2]);  # get result
      agena := f(argu!0);
      if maple = 'singularityencountered' or
         maple = 'divisionbyzero' then maple := undefined fi;
      if maple = 'true' then maple := true fi;
      if maple = 'false' then maple := false fi;
      if type(maple) <> boolean then
         if isdiff(agena, maple) then
            print('>>> ', file, argu, agena, maple)
         fi
      else
         if agena <> maple then
            print('>>> ', file, argu, agena, maple)
         fi
     fi
   od;
end;


cmplxtest2 := proc(file::string, f::procedure) is
   local s, se, argu, maple, agena;
   for keys j in io.lines(file) do
      if j[1] = '#' then skip fi;
      s := replace(j, ' ', '');  # delete white spaces
      se := s split ';';         # split into its components
      argu := tonumber(se[1]);   # get argument
      maple := tonumber(se[2]);  # get result
      agena := f(argu);
      if maple = 'singularityencountered' or
         maple = 'divisionbyzero' then maple := undefined fi;
      if maple = 'true' then maple := true fi;
      if maple = 'false' then maple := false fi;
      if type(maple) <> boolean then
         if isdiff(agena, maple) then
            print('>>> ', file, argu, agena, maple)
         fi
      else
         if agena <> maple then
            print('>>> ', file, argu, agena, maple)
         fi
     fi
   od;
end;


cmplxtest('math3.tst', << x -> x*exp(x)*sin(x)-ln(cos(x)) >>);
cmplxtest('math4.tst', << x -> x*exp(x)*ln(sin(x))-ln(cos(x)) >>);
cmplxtest('math5.tst', << x -> x*exp(x)*ln(tan(x))-ln(cos(x)) >>);
cmplxtest('math6.tst', << x -> x*exp(x)*ln(sinh(x))-ln(cos(x)) >>);
cmplxtest('math7.tst', << x -> x*exp(x)*ln(cosh(x))-ln(cos(x)) >>);
cmplxtest('math8.tst', << x -> x*exp(x)*ln(tanh(x))-ln(cos(x)) >>);
cmplxtest('math9.tst', << x -> x*exp(x)*ln(arcsin(x))-ln(cos(x)) >>);
cmplxtest('matha.tst', << x -> x*exp(x)*ln(arccos(x))-ln(cos(x)) >>);
cmplxtest('mathb.tst', << x -> x*exp(x)*ln(arctan(x))-ln(cos(x)) >>);
cmplxtest('mathc.tst', << x -> x*exp(x)*ln(sqrt(x))-ln(cos(x)) >>);
cmplxtest('mathd.tst', << x -> x*exp(x)*ln(lngamma(x))-ln(cos(x)) >>);
cmplxtest2('mathe.tst', << x -> calc.xpdiff(<< x -> sin(x) >>, x) >>);
cmplxtest2('mathe.tst', << x -> calc.diff(<< x -> sin(x) >>, x) >>);
cmplxtest('mathf.tst', << x -> (x+I*(x-1) + (x+1)+I*(x-1)) >>);
cmplxtest('mathg.tst', << x -> (x+I*(x-1) - (x+1)+I*(x-1)) >>);
cmplxtest('mathh.tst', << x -> ((x+I*(x-1)) * ((x+3)+I*(x-2))) >>);
cmplxtest('mathi.tst', << x -> ((x+I*(x-1)) / ((x+3)+I*(x-2))) >>);
cmplxtest2('mathj.tst', << x -> ((x+I*(x-1)) ** (round(x) % 32)) >>);

print("\nChecking prime functions\n");

try  # prime bug in Agena 10.5 and earlier
   if math.isprime(4486052565) <> false or math.isprime(76234567895) <> false then
      print('error in math.isprime, 1')
   fi;
   if math.prevprime(573259391) <> 573259369 then
      print('error in math.isprime, 2')
   fi;
   if math.nextprime(573259369) <> 573259391 then
      print('error in math.isprime, 3')
   fi;
catch errmsg then
   print('error in math.isprime, 4:', errmsg)
yrt;

L := {};
r := 0;
f := binio.open('isprime.bin');
a := binio.readchar(f);
while a <> null do
   inc r, a;
   insert r into L;
   a := binio.readchar(f)
od;
binio.close(f);

for i from -2 to 1m do
   temp := math.isprime(i);
   if temp and not(i in L) then
      print('error in math.isprime, 4(' & i & ')')
   fi
od;
clear L, r, f, temp;

print("\nChecking linalg package\n");

if unassigned(linalg) then
   import linalg;
fi;

A := linalg.matrix([1, 1, 1], [1, 2, 3], [2, 3, 6])
b := linalg.vector(-6, -10, -18)
if linalg.gsolve(A, b) <> linalg.vector(-3, -2, -1) then
   print('error in linalg.gsolve, 1')
fi;

v := linalg.vector(1.1, 2.2, 3.3);

if v + v <> linalg.vector(2.2, 4.4, 6.6) then
   print('Error in `linalg.add')
fi;

if 2*v - v <> v then
   print('Error in `linalg.scalarmul/linalg.sub')
fi;

if linalg.scalarmul(v, 2) <> linalg.scalarmul(2, v) then
   print('Error in `linalg.scalarmul, 1')
fi;

a := linalg.matrix([0, 1, 0], [2, 3, 1]);

b := linalg.matrix([2, 0], [1, 2], [1, 1]);

M := linalg.matrix([1, 1, 1, 0], [2, 1, 1, 1], [1, 2, 1, 15]);

A := linalg.matrix([1, 1, 1], [2, 1, 1], [1, 2, 1]);

v := linalg.vector(0, 1, 15);

m := linalg.matrix([1, 2, 0], [-1, 1, 1], [1, 2, 3]);

to 10 do
   if not(a * b == linalg.matrix([1, 2], [8, 7])) then
      print('Error in linalg.mmul, 1')
   fi;
   if not(b * a == linalg.matrix([0, 2, 0], [4, 7, 2], [2, 4, 1])) then
      print('Error in linalg.mmul, 2')
   fi;
   if not(linalg.identity(3) == linalg.matrix([1, 0, 0], [0, 1, 0], [0, 0, 1])) then
      print('Error in linalg.identity or linalg.meq, 1')
   fi;
   if not(linalg.identity(3) * linalg.matrix([1, 0, 0], [0, 1, 0], [0, 0, 1]) == linalg.identity(3)) then
      print('Error in linalg.identity or linalg.meq, 2')
   fi;
   if not(2 * linalg.identity(2) == linalg.matrix([2, 0], [0, 2])) then
     print('Error in `linalg.scalarmul/linalg.meq`')
   fi;
   if linalg.det(linalg.matrix([1, 2, 3], [4, 5, 6], [7, 8, 9])) <> 0 then
     print('Error in `linalg.det`, 1')
   fi;
   if linalg.det(linalg.matrix([1, 2, 3], [4, 5, 6], [7, 8, 10])) <> -3 then
      print('Error in `linalg.det`, 2')
   fi;
   if linalg.det(linalg.matrix([1, 2, -1, 10, 11], [4, 5, 6, 12, 13],
      [7, 8, 9, 14, 15], [2, 3, 4, 5, 6], [10, 11, 12, 13, 16])) <> -136 then
      print('Error in `linalg.det`, 3')
   fi;
   if not(linalg.inverse(m) ~=
      linalg.matrix([0.11111111111111, -0.66666666666667, 0.22222222222222],
                    [0.44444444444444,  0.33333333333333, -0.11111111111111],
                    [-0.33333333333333, 0, 0.33333333333333])) then
      print('Error in `linalg.inverse`')
   fi;
   if not(linalg.transpose(m) ==
      linalg.matrix([1, -1, 1],
                    [2,  1, 2],
                    [0,  1, 3])) then
      print('Error in `linalg.transpose, 1`')
   fi;
   if not(linalg.transpose(linalg.matrix([1, 2, 0], [-1, 1, 1], [1, 2, 3], [0, 0, 0], [8, 8, 8])) ==
      linalg.matrix([1, -1, 1, 0, 8], [2,  1, 2, 0, 8], [0,  1, 3, 0, 8] )) then
      print('Error in `linalg.transpose, 2`')
   fi;
   if not(linalg.diagonal(linalg.vector(1, 1, 3)) == linalg.matrix([1, 0, 0], [0, 1, 0], [0, 0, 3])) then
      print('Error in `linalg.diagonal`')
   fi;
   if not(linalg.getdiagonal(m) == linalg.vector(1, 1, 3)) then
      print('Error in `linalg.getdiagonal`')
   fi;
   if linalg.trace(m) <> 5 then
      print('Error in `linalg.trace`')
   fi;
   if not(linalg.gsolve(M) == linalg.gsolve(A, v)) then
      print('Error in `linalg.gsolve`, 1')
   fi;
   if linalg.gsolve(linalg.matrix([1, 1, 1], [2, 2, 5], [4, 4, 8]), linalg.vector(-1, -8, -14)) <> fail then  # either no or infinite solutions
      print('Error in `linalg.gsolve, 2`')
   fi;
   if linalg.gsolve(linalg.matrix([1, -2, -3], [3, 1, 1], [11, -1, -3]), linalg.vector(-1, 4, 10)) <> infinity then
      print('Error in `linalg.gsolve, 3`')
   fi;
   if linalg.gsolve(linalg.matrix([2, 4, 8], [1, 1, 1], [0, 2, 6]), linalg.vector(8, 3, 2)) <> infinity then
      print('Error in `linalg.gsolve, 4`')
   fi;
   if linalg.gsolve(linalg.matrix([2, 4, 8], [1, 1, 1], [0, 2, 6]), linalg.vector(8, 3, 1)) <> undefined then
      print('Error in `linalg.gsolve, 5`')
   fi;
   if not(linalg.gsolve(linalg.matrix([-1, 1, 2], [3, -1, 1], [-1, 3, 4]), linalg.vector(1, 1, 1)) ~= linalg.vector(-0.2, -0.8, 0.8)) then
      print('Error in `linalg.gsolve, 6`')
   fi;
   if not(linalg.gsolve(linalg.matrix([3, 1, -1], [0, 2, 1], [0, 0, 3]), linalg.vector(0, 12, 6)) ~= linalg.vector(-1, 5, 2)) then
      print('Error in `linalg.gsolve, 7`')
   fi;
   if not(linalg.gsolve(linalg.matrix([2, 4, 8], [1, 1, 1], [1, 3, 1]), linalg.vector(14, 3, 5)) ~= linalg.vector(1, 1, 1)) then
      print('Error in `linalg.gsolve, 8`')
   fi
   xx, xy, xz:= linalg.ludecomp(linalg.matrix([1, 2, 3], [-1, 2, 1], [0, 1, 2]), 3);
   if not(xx ~= linalg.matrix([1, 2, 3], [-1, 4, 4], [0, 0.25, 1])) or
      not(xy ~= linalg.vector([1, 2, 3])) or xz <> 1 then
      print('Error in `linalg.ludecomp`')
   fi;
   if linalg.submatrix(linalg.matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]), 2:3, 1:3) <> linalg.matrix([2, 3], [5, 6], [8, 9]) then
      print('Error in `linalg.submatrix, 1`')
   fi;
   if linalg.submatrix(linalg.matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]), 2:2, 1:3) <> linalg.vector([2, 5, 8]) then
      print('Error in `linalg.submatrix, 2`')
   fi;
   xx, xy := linalg.rref(linalg.matrix([2, 4, 8, 8], [1, 1, 1, 3], [0, 2, 6, 1]));
   if not(xx ~= linalg.matrix([1, 0, -2, 0], [0, 1,  3, 0], [0, 0, 0, 1])) or xy <> fail then
      print('Error in `linalg.rref, 1`')
   fi;
   xx, xy := linalg.rref(linalg.matrix([2, 4, 8, 8], [1, 1, 1, 3], [0, 0, 0, 0]));
   if not(xx ~= linalg.matrix([1, 0, -2, 2], [0, 1,  3, 1], [0, 0, 0, 0])) or xy <> fail then
      print('Error in `linalg.rref, 2`')
   fi;
   xx, xy := linalg.rref(linalg.matrix([3, 1, -1, 0], [0, 2, 1, 12], [0, 0, 3, 6]));
   if not(xx ~= linalg.matrix([1, 0, 0, -1], [0, 1, 0,  5], [0, 0, 1, 2])) or xy <> null then
      print('Error in `linalg.rref, 3`')
   fi;
   xx, xy := linalg.rref(linalg.matrix([-1, 1, 2, 1], [3, -1, 1, 1], [-1, 3, 4, 1]));
   if not(xx ~= linalg.matrix([1, 0, 0, -0.2], [0, 1, 0, -0.8], [0, 0, 1, 0.8])) or xy <> null then
      print('Error in `linalg.rref, 4`')
   fi;
   xx, xy := linalg.rref(linalg.matrix([1,-2,-3,-1], [3, 1, 1, 4], [11, -1, -3, 10]));
   if not(xx ~= linalg.matrix([1, 0, -0.14285714285714, 1], [0, 1, 1.4285714285714, 1], [0, 0, 0, 0])) or xy <> fail then
      print('Error in `linalg.rref, 5`')
   fi;
   xx, xy := linalg.rref(linalg.matrix([1, 2, 2], [0, 1, 0], [0, 0, 1]), linalg.vector(9, 1, 2));
   if not(xx ~= linalg.matrix([1, 0, 0, 3], [0, 1, 0, 1], [0, 0, 1, 2])) or xy <> null then
      print('Error in `linalg.rref, 6`')
   fi;
   if not(linalg.backsub(linalg.matrix([1, 2, 2], [0, 1, 0], [0, 0, 1]), linalg.vector(9, 1, 2)) ~= linalg.vector(3, 1, 2)) then
      print('Error in `linalg.backsub, 1`')
   fi;
   if not(linalg.forsub(linalg.matrix([1, 0, 0], [0, 1, 0], [2, 2, 1]), linalg.vector(9, 1, 2)) ~= linalg.vector(9, 1, -18)) then
      print('Error in `linalg.forsub, 1`')
   fi
od;

print("\nChecking stats package\n");

# stats package
A1 := [1., 1.414213562, 1.732050808, 2., 2.236067978, 2.449489743, 2.645751311];
A2 := [2.449489743, 2.645751311];
A3 := [2.449489743];
A5 := totable(nseq(<< x -> 0.1 >>, 1, 10));
A6 := [undefined, 1., 1.414213562, undefined, 1.732050808, undefined, 2., undefined, 2.236067978, undefined, 2.449489743, undefined, 2.645751311, undefined];
A7 := totable(nseq(<< x -> sin(x) >>, 1, 100));

lew := read('lew_nist.agb');

stattest := proc(A1, A2, A3, A5, A6, A7) is
   if not approx(stats.sd(A1), 0.5412573272) then print('error in stats.sd, 1') fi;
   if not approx(stats.sd(A2), .09813078401) then print('error in stats.sd, 2') fi;
   if not approx(stats.sd(A7), 0.70900088711577) then print('error in stats.sd, 3') fi;
   if not approx(stats.sd(A7, true), 0.71257270260368) then print('error in stats.sd, 4') fi;
   if stats.sd([]) <> fail then print('error in stats.sd, 5') fi;

   if not approx(stats.ad(A1), 0.4656681477) then print('error in stats.ad, 1') fi;
   if not approx(stats.ad(A2), .9813078400e-1) then print('error in stats.ad, 2') fi;
   if stats.ad(A3) <> fail then print('error in stats.ad, 3') fi;
   if stats.ad([]) <> fail then print('error in stats.ad, 4') fi;

   if not approx(stats.gmean(A1), 1.838481754) then print('error in stats.gmean, 1') fi;
   if not approx(stats.gmean(A2), 2.545729895) then print('error in stats.gmean, 2') fi;
   if stats.gmean([]) <> fail then print('error in stats.gmean, 3') fi;

   if not approx(stats.hmean(A1), 1.742210832663344) then print('error in stats.hmean, 1') fi;
   if not approx(stats.hmean(A2), 2.543840666276475) then print('error in stats.hmean, 2') fi;
   if stats.hmean([]) <> fail then print('error in stats.hmean, 3') fi;

   if not approx(stats.mean(A1), 1.925367629) then print('error in stats.mean, 1') fi;
   if not approx(stats.mean(A2), 2.547620527) then print('error in stats.mean, 2') fi;
   if stats.mean([]) <> fail then print('error in stats.mean, 3') fi;

   if not approx(stats.amean(A1), 1.925367629) then print('error in stats.amean, 1') fi;
   if not approx(stats.amean(A2), 2.547620527) then print('error in stats.amean, 2') fi;
   if stats.amean([]) <> fail then print('error in stats.amean, 3') fi;

   if not approx(stats.qmean(A1), 2.000000000) then print('error in stats.qmean, 1') fi;
   if not approx(stats.qmean(A2), 2.549509757) then print('error in stats.qmean, 2') fi;
   if stats.qmean([]) <> fail then print('error in stats.qmean, 3') fi;

   if not approx(stats.median(A1), 2.000000000) then print('error in stats.median, 1') fi;
   if not approx(stats.median(A2), 2.547620527) then print('error in stats.median, 2') fi;
   if stats.median([]) <> fail then print('error in stats.median, 3') fi;

   if not approx(stats.var(A1), 0.2929594942) then print('error in stats.var, 1') fi;
   if not approx(stats.var(A2), .9629650770e-2) then print('error in stats.var, 2') fi;
   if not approx(stats.var(A7), 0.50268225793095) then print('error in stats.var, 3') fi;
   if not approx(stats.var(A7, true), 0.50775985649591) then print('error in stats.var, 4') fi;
   if stats.var([]) <> fail then print('error in stats.var, 5') fi;

   if not approx(stats.skewness(A1), -.3468060723) then print('error in stats.skewness, 1') fi;
   if not approx(stats.skewness(A2), 0) then print('error in stats.skewness, 2') fi;
   if stats.skewness(A3) <> fail then print('error in stats.skewness, 3') fi;
   if stats.skewness([]) <> fail then print('error in stats.skewness, 4') fi;

   if not approx(stats.moment(A1, 5, 1), 3.225641950) then print('error in stats.moment, 1') fi;
   if not approx(stats.moment(A2, 5, 1), 9.235813300) then print('error in stats.moment, 2') fi;
   if stats.moment([], 5, 1) <> fail then print('error in stats.moment, 3') fi;

   if not approx(stats.iqmean(A1), 1.845583087) then print('error in stats.iqmean, 1') fi;
   if not approx(stats.iqmean(A2), 2.449489743) then print('error in stats.iqmean, 2') fi;
   if stats.iqmean(A3) <> fail then print('error in stats.iqmean, 3') fi;

   local B1, B2;
   if type(A1) = table then
      B1, B2 := [1, 2.645751311], [2.449489743, 2.645751311]
   else
      B1, B2 := seq(1, 2.645751311), seq(2.449489743, 2.645751311)
   fi;
   if stats.minmax(A1) <> B1 then print('error in stats.minmax, 1') fi;
   if stats.minmax(A2) <> B2 then print('error in stats.minmax, 2') fi;
   if stats.minmax(A3) <> fail then print('error in stats.minmax, 3') fi;
   if stats.minmax([]) <> fail then print('error in stats.minmax, 4') fi;

   local A4 := [6, 12, 12, 12, 10, 5];
   if type(A1) = table then
      B1 :=[12]
   else
      A4 := toseq(A4);
      B1 := seq(12)
   fi;
   if stats.mode(A4) <> B1 then print('error in stats.mode/1') fi;
   if toset(stats.mode(A1)) <> {2.645751311, 1.732050808, 1, 2, 2.449489743, 1.414213562, 2.236067978} then print('error in stats.mode/2') fi;

   local r, s, s1, s2, s3, s4, s5;
   r := stats.cumsum(A5);
   s := seq(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1);
   s1 := seq(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, undefined, 1, undefined);
   local flag := false;
   for i to size s do
      if abs(r[i] - s[i]) > Eps then
         flag := true; print('error in stats.cumsum, 1 ' & r[i] & ' ' & s[i] & ' ' & abs(r[i] - s[i]))
      fi;
   od;
   if flag then
      print('error in stats.cumsum, 1')
   fi;
   if stats.cumsum([]) <> fail or stats.cumsum(seq()) <> fail then
      print('error in stats.cumsum, 2')
   fi;
   s2 := seq(0.1, 0.3, 0.6, 1, 1.5, 2.1, 2.8, 3.6, 4.5, undefined, 5.5, undefined);
   s3 := stats.cumsum(s1);
   for i to size s2 do
      if approx(s2[i], s3[i]) <> true then print('error in stats.cumsum, 3') fi
   od;
   s2 := totable(s2);
   s3 := stats.cumsum(totable(s1));
   for i to size s2 do
      if approx(s2[i], s3[i]) <> true then print('error in stats.cumsum, 4') fi
   od;
   if stats.sum(s) <> 5.5 then
      print('error in stats.sum, 1')
   fi;
   if stats.sum(totable(s)) <> 5.5 then
      print('error in stats.sum, 2')
   fi;
   s4 := [4, 3, 2, 1];
   if stats.sorted(s4) <> sorted(s4) then
      print('error in stats.sorted, 1')
   fi;
   s4 := toseq(s4);
   if stats.sorted(s4) <> sorted(s4) then
      print('error in stats.sorted, 2')
   fi;
   s5 := [4, 3, 2, 1, -1, -2, -4, -5, undefined];
   if [stats.colnorm(s5)] <> [5, -5] then
      print('error in stats.colnorm, 1')
   fi;
   s5 := toseq(s5);
   if [stats.colnorm(s5)] <> [5, -5] then
      print('error in stats.colnorm, 2')
   fi;
   s5 := totable(s5);
   if stats.rownorm(s5) <> 22 then
      print('error in stats.rownorm, 1')
   fi;
   s5 := toseq(s5);
   if stats.rownorm(s5) <> 22 then
      print('error in stats.rownorm, 2')
   fi;
   s5 := totable(s5);
   if stats.issorted(s5) then
      print('error in stats.issorted, 1')
   fi;
   s5 := toseq(s5);
   if stats.issorted(s5) then
      print('error in stats.issorted, 2')
   fi;
   local s6 := seq(4, 3, 2, 1, 1, 2, 4, 5)
   if stats.countentries(s6) <> [2, 2, 1, 2, 1] then
      print('error in stats.countentries, 1')
   fi;
   s6 := totable(s6);
   if stats.countentries(s6) <> [2, 2, 1, 2, 1] then
      print('error in stats.countentries, 2')
   fi;
   if stats.smallest(s6, 3) <> 2 then
      print('error in stats.smallest')
   fi;
   if approx(stats.pdf(36, 35, 2), 0.17603266338215) <> true then
      print('error in stats.pdf')
   fi;
   if approx(stats.ndf(2) * stats.nde(36, 35, 2), 0.17603266338215) <> true then
      print('error in stats.ndf or stats.nde')
   fi;
   if approx(stats.cdf(1, 3, 2, 6), 0.13236766522181) <> true then
      print('error in stats.cdf, 1')
   fi;
   if approx(stats.cdf(-1, 3, 2, 6), 0.25764629388492) <> true then
      print('error in stats.cdf, 2')
   fi;
   if approx(stats.sma(s, 4, 3, true), 0.4) <> true then
      print('error in stats.sma, 1')
   fi;
   if protect(stats.sma, s1, 10, 3, true) <> lasterror then
      print('error in stats.sma, 2')
   fi;
   if approx(stats.smm(s, 4, 3, true), 0.4) <> true then
      print('error in stats.smm, 1')
   fi;
   if protect(stats.smm, s1, 10, 3, true) <> lasterror then
      print('error in stats.smm, 2')
   fi;
   f := stats.gsma(s, 4, 3, true);
   for i from 4 to 7 do
      if f() <> stats.sma(s, i, 3, true) then
         print('Error in stats.gsma/1, ' & i-3)
      fi
   od;
   f := stats.gsmm(s, 4, 3, true);
   for i from 4 to 7 do
      if f() <> stats.smm(s, i, 3, true) then
         print('Error in stats.gsmm/1, ' & i-3)
      fi
   od;
   f := stats.gsma(s, 4, 3);
   for i from 4 to 7 do
      if f() <> stats.sma(s, i, 3) then
         print('Error in stats.gsma/2, ' & i-3)
      fi
   od;
   f := stats.gsmm(s, 4, 3);
   for i from 4 to 7 do
      if f() <> stats.smm(s, i, 3) then
         print('Error in stats.gsmm/3, ' & i-3)
      fi
   od;
   if stats.trimmean(seq(11.4, 17.3, 21.3, 25.9, 40.1), 0.4) <> 21.5 then
      print('Error in stats.trimmean')
   fi;
   if not approx(stats.herfindahl(seq(2300, 2500, 1900, 1800, 3400, 2700, 2900, 2100)), 0.0060688998631522) then
      print('Error in stats.herfindahl')
   fi;
   s := seq(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1);
   if not(sqrt(stats.fsum(<< x, y -> x**2 + y >>, s, 1, size s, 0)/size s - stats.amean(s)**2) ~= stats.sd(s)) then
      print('Error in stats.fsum')
   fi;
   if not(stats.fprod(<< x, y -> x**2 + y >>, s, 1, size s, 0) ~= 1.316818944e-007) then
      print('Error in stats.fprod')
   fi;
   if not(stats.acf(lew, 0) ~= 1) then
      print('Error in stats.acf')
   fi;
   if not(stats.acv(lew, 0) ~= seq(1)) then
      print('Error in stats.acv')
   fi;
   if not(stats.neighbours(seq(seq(1, 1, 1), seq(1, 1, 2)), 1, 2, 1) ~= seq(seq(1, 1, 2))) then
      print('Error in stats.neighbours, 1')
   fi;
   if not(stats.neighbours(seq(linalg.vector(1, 1, 1), linalg.vector(1, 1, 2)), 1, 2, 1) ~= seq(linalg.vector(1, 1, 2))) then
      print('Error in stats.neighbours, 2')
   fi;
   a, b := stats.meanvar(A7);
   if a ~<> -0.0012717101366042 or b ~<> 0.50268225793095 then
      print('Error in stats.meanvar, 1')
   fi;
   a, b := stats.meanvar(A7, true);
   if a ~<> -0.0012717101366042 or b ~<> 0.50775985649591 then
      print('Error in stats.meanvar, 2')
   fi;
   if stats.covar(A7, nseq( << x -> cos x >>, 1, 100)) ~<> -0.0013668298931265 then
      print('Error in stats.covar')
   fi;
   if stats.standardise(totable(nseq( << x -> sin x >>, 1, 5))) ~<> [1.0338280579828, 1.120800914479, 0.1357782239415, -1.0156146946144, -1.2747925017889] then
      print('Error in stats.standardise')
   fi
end;


if stats.chauvenet(seq()) <> fail or stats.chauvenet(seq(1)) <> fail then
   print('Error in stats.chauvenet, 1')
fi;


s := seq(8.02, 8.16, 3.97, 8.64, 0.84, 4.46, 0.81, 7.74, 8.78, 9.26, 20.46, 29.87, 10.38, 25.71)

if stats.chauvenet(s) <> seq(29.87, 25.71, 20.46) then
   print('Error in stats.chauvenet, 2')
fi;

if stats.chauvenet(s, 29.87) <> true or stats.chauvenet(s, 25.71) <> true or stats.chauvenet(s, 20.46) <> false then
   print('Error in stats.chauvenet, 3')
fi;


s := seq();

s1 := seq();

for i to 100 do
   insert i into s
od;

for i from 100 downto 1 do
   insert i into s
od;

for i to 100 do
   to 2 do insert i into s1 od
od

if stats.sorted(s) <> s1 then
   print('Error in stats.sorted, 1')
fi;

if stats.sorted(s, 'heapsort') <> s1 then
   print('Error in stats.sorted, 2')
fi;

if stats.sorted(s, 'quicksort') <> s1 then
   print('Error in stats.sorted, 3')
fi;

if stats.sorted(s, 'nrquicksort') <> s1 then
   print('Error in stats.sorted, 4')
fi;

if stats.sorted(s, 'introsort') <> s1 then
   print('Error in stats.sorted, 5')
fi;

if stats.gini(s) <> 0.33 then
   print('Error in stats.gini, 1')
fi;

if stats.gini(s1, 'sorted') <> 0.33 then
   print('Error in stats.gini, 2')
fi;

if not(stats.chisquare(2, 3) ~= .20755374871029735167) then
   print('Error in stats.chisquare, 1')
fi;

if not(stats.normald(2, 3, 2) ~= .17603266338214973889) then
   print('Error in stats.normald, 1')
fi;

if not(stats.fratio(2, 3, 2) ~= .12178482240718668470) then
   print('Error in stats.fratio, 1')
fi;

if not(stats.studentst(-2, 3) ~= .67509660663892904023e-1) then
   print('Error in stats.studentst, 1')
fi;

if not(stats.kurtosis([0,1,1,10,0,1]) ~= 4.0872862733815804485) then
   print('Error in stats.kurtosis, 1')
fi;

if not(stats.kurtosis([0,1,1,10,11,0,1]) ~= 1.9281632088969823101) then
   print('Error in stats.kurtosis, 2')
fi;


stattest(A1, A2, A3, A5, A6, A7);
stattest(toseq(A1), toseq(A2), toseq(A3), toseq(A5), toseq(A6), toseq(A7));


for i from 10 to 50 by 10 do
   a := stats.prange(seq(1, 2), i);
   a := stats.prange(seq(1, 2, 3), i);
   a := stats.prange(seq(1, 2, 3, 4), i);
   b := stats.winsor(seq(1, 2), i);
   b := stats.winsor(seq(1, 2, 3), i);
   b := stats.winsor(seq(1, 2, 3, 4), i);
od;


create seq s(100);

for i to 30 do insert 1 into s od;

for i to 60 do insert 2 into s od;

insert 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 into s;

to 10 do
   if stats.weights(1, 2, 3, 4, 5, 6, 1:30, 2:60, 7, 8, 9, 10) <> s then
      error('Error in stats.weights')
   fi
od;

clear a, b;

print("\nChecking calc package\n");

to 10 do
   f := calc.linterp([ 0:0, 1:3, 2:1, 3:3 ]);
   if f(1) <> 3 then
      print('Error in calc.linterp, 1')
   fi;

   f := calc.linterp(seq(0:0, 1:3, 2:1, 3:3));
   if f(1) <> 3 then
      print('Error in calc.linterp, 2')
   fi;

   f := calc.neville([ 0:0, 1:3, 2:1, 3:3 ]);
   if f(1) <> 3 then
      print('Error in calc.neville, 1')
   fi;

   f := calc.neville(seq(0:0, 1:3, 2:1, 3:3));
   if f(1) <> 3 then
      print('Error in calc.neville, 2')
   fi;

   f := calc.linterp([1:0.1, 2:3, 3:(-4), 4:1, 5.1:0]);
   if not(approx(f(1), 0.1) or
      approx(f(2), 3) or
      approx(f(3), -4) or
      approx(f(4), 1) or
      approx(f(5), 1.8911512562968) or
      approx(f(5.1), 0) or
      approx(f(6), -55.444243718516)) then
         print('Error in calc.linterp, 3')
   fi;

   f := calc.linterp(seq(1:0.1, 2:3, 3:(-4), 4:1, 5.1:0));
   if not(approx(f(1), 0.1) or
      approx(f(2), 3) or
      approx(f(3), -4) or
      approx(f(4), 1) or
      approx(f(5), 1.8911512562968) or
      approx(f(5.1), 0) or
      approx(f(6), -55.444243718516)) then
         print('Error in calc.linterp, 4')
   fi;

   f := calc.interp([ 0:0, 1:3, 2:1, 3:3 ]);
   if f(1) <> 3 then
      print('Error in calc.interp')
   fi;

   f := calc.polyfit([1:0.1, 2:3, 3:(-4), 4:1, 5.1:0], 3);
   if not(approx(f[1], 0.26689930948382) or
      approx(f[2], -2.1617356885251) or
      approx(f[3], 4.5685918608638) or
      approx(f[4], -1.9802624870949)) then
         print('Error in calc.polyfit')
   fi;

   f := calc.polyfit(seq(1:0.1, 2:3, 3:(-4), 4:1, 5.1:0), 3);
   if not(approx(f[1], 0.26689930948382) or
      approx(f[2], -2.1617356885251) or
      approx(f[3], 4.5685918608638) or
      approx(f[4], -1.9802624870949)) then
         print('Error in calc.polyfit')
   fi;

   if not(calc.gtrap(<< x -> sin(x) >>, 0, Pi) ~= 1.9999999938721) then
      print('Error in calc.gtrap')
   fi;

   if not(calc.integ(<< x -> sin(x) >>, 0, Pi) ~= 2) then
      print('Error in calc.integ')
   fi;

   if not(calc.simaptive(<< x -> sin(x) >>, 0, Pi) ~= 2.0000000003246) then
      print('Error in calc.simaptive')
   fi;

   a, b := calc.minimum(<< x -> sin x >>, -10, 10);
   if a ~<> seq(-7.8539815784491, -1.5707963196727, 4.7123889665853) or
      b ~<> seq(-7.8539815784491) then
      print('Error in calc.minimum')
   fi;
   clear a, b
od;

create table f;
for i to 3 do f[i] := [] od;

f[1, 0] := << x -> x^3 - x^2 + x - 1 >>;
f[1, 1] := << x -> 3*x^2 - 2*x + 1 >>;
f[1, 2] := << x -> 6*x - 2 >>;
f[1, 3] := << x -> 6 >>;
f[2, 0] := << x -> sin x >>;
f[2, 1] := << x -> cos x >>;
f[2, 2] := << x -> -sin x >>;
f[2, 3] := << x -> -cos x >>;
f[3, 0] :=  << x -> ln x >>;
f[3, 1] := << x -> 1/x >>;
f[3, 2] := << x -> -1/(x*x) >>;
f[3, 3] := << x -> 2/(x*x*x) >>;

for k to 3 do
   for i from 0 to 3 do
      #c := calc.chebycoeffs(f[k, 0], 1, 10, 50);
      g := calc.cheby(f[k, 0], 1, 10, 50, deriv = i);
      for x from 2 to 9 do
         if f[k, i](x) ~<> g(x) then
            print('error in cheby/cheby, 1', i, x, f[k, i](x) - g(x))
         fi
      od;
   od;
od;

if calc.gaussian(2, 5, 10, 15) ~<> 4.3371423659171 then
   print('error in calc.gaussian, 1')
fi;

if calc.gaussian(2!1, 5, 10, 15) ~<> 4.3440438651045+0.15452001307634*I then
   print('error in calc.gaussian, 2')
fi;

if calc.eulerdiff( << x -> ln x >>, 3) ~<> 1/3 then
   print('error in calc.eulerdiff, 1')
fi;

clear c, f, g, i, k, x;


print("\nChecking cordic package\n");

for i to 10 do
   for x from -100 to 100 by .1 do
      if not(ln(x) ~= cordic.cln(x)) then
         print('Error in `cordic.cln`')
      fi;
      if not(exp(x) ~= cordic.cexp(x)) then
         print('Error in `cordic.cexp`')
      fi;
      if not(sqrt(x) ~= cordic.csqrt(x)) then
         print('Error in `cordic.csqrt`')
      fi;
      if not(cbrt(x) ~= cordic.ccbrt(x)) then
         print('Error in `cordic.ccbrt`')
      fi;
      if not(sin(x) ~= cordic.csin(x)) then
         print('Error in `cordic.csin`')
      fi;
      if not(cos(x) ~= cordic.ccos(x)) then
         print('Error in `cordic.ccos`')
      fi;
      if not(tan(x) ~= cordic.ctan(x)) then
         print('Error in `cordic.ctan`')
      fi;
      if not(sinh(x) ~= cordic.csinh(x)) then
         print('Error in `cordic.csinh`')
      fi;
      if not(cosh(x) ~= cordic.ccosh(x)) then
         print('Error in `cordic.ccosh`')
      fi;
      if not(tanh(x) ~= cordic.ctanh(x)) then
         print('Error in `cordic.ctanh`')
      fi;
      if not(arcsin(x) ~= cordic.carcsin(x)) then
         print('Error in `cordic.carcsin`')
      fi;
      if not(arccos(x) ~= cordic.carccos(x)) then
         print('Error in `cordic.carccos`', x, arccos(x), cordic.carccos(x))
      fi;
      if not(arctan2(x, x+1) ~= cordic.carctan2(x, x+1)) then
         print('Error in `cordic.carctan2`')
      fi;
      if not(hypot(x, x+1) ~= cordic.chypot(x, x+1)) then
         print('Error in `cordic.chypot`')
      fi;
      if not(arctanh(x) ~= cordic.carctanh(x)) then
         print('Error in `cordic.carctanh`')
      fi;
   od
od;


#print('before polygen');

for i to 10 do
   g := calc.polygen(5, 4, 3, 2, 1);
   if g(0) <> 1 or g(1) <> 15 or g(2) <> 129 then
      print('Error in calc.polygen, 1')
   fi
od;

to 10 do
   s := nseq(1, 256);
   reverse(s);
   g := calc.polygen(unpack(s));
   if g(0) ~<> 1 or g(0.5) ~<> 4 or g(0.1) ~<> 1.2345679012346 or g(0.2) ~<> 1.5625 then
      print('Error in calc.polygen, 2')  # checked with Maple V Release 4
   fi
od;


print("\nChecking divs package\n");

if divs.divs(1, 2, 3) + divs.divs(1, 3) <> divs.divs(2) then
   print('Error in divs/+, 1')
fi;

if divs.divs(1, 2) + divs.divs(1, 3) <> divs.divs(5, 6) then
   print('Error in divs/+, 2')
fi;

if divs.divs(1, 2) * divs.divs(1, 3) <> divs.divs(1, 6) then
   print('Error in divs/*, 1')
fi;

if 2 * divs.divs(1, 3) <> divs.divs(2, 3) then
   print('Error in divs/*, 2')
fi;

if divs.divs(1, 3) / divs.divs(2, 3) <> divs.divs(1, 2) then
   print('Error in divs//')
fi;

if sign divs.divs(1, -3) <> -1 then
   print('Error in divs/sign')
fi;

if sqrt divs.divs(2, 3) <> divs.divs(175937443, 215478481) then
   print('Error in divs/sqrt')
fi;

if not(ln divs.divs(2, 3) ~= divs.divs(-31683730, 78141693)) then
   print('Error in divs/ln')
fi;

if exp divs.divs(2, 3) <> divs.divs(215873355, 110833076) then
   print('Error in divs/exp')
fi;

if sin divs.divs(2, 3) <> divs.divs(70029431, 113248465) then
   print('Error in divs/sin')
fi;

if cos divs.divs(2, 3) <> divs.divs(167593363, 213253696) then
   print('Error in divs/cos')
fi;

if tan divs.divs(2, 3) <> divs.divs(223635209, 284218377) then
   print('Error in divs/tan')
fi;

if arctan divs.divs(2, 3) <> divs.divs(97020125, 164999482) then
   print('Error in divs/arctan')
fi;

if divs.divs(2, 4) < divs.divs(2, 3) <> true then
   print('Error in divs/relation, 1')
fi;

if divs.divs(2, 4) > divs.divs(2, 3) <> false then
   print('Error in divs/relation, 2')
fi;

if divs.divs(2, 4) = divs.divs(1, 2) <> true then
   print('Error in divs/relation, 3')
fi;

if divs.todec(divs.divs(1, 2)) <> 0.5 then
   print('Error in divs.todec')
fi;

a, b := divs.todiv(0.5);

if a <> divs.divs(1, 2) or b <> 0 then
   print('Error in divs.todiv')
fi;

if abs(divs.divs(-1, 3)) <> divs.divs(1, 3) then
   print('Error in divs/abs')
fi;

if divs.denom(divs.divs(-1, 3)) <> 3 then
   print('Error in divs.denom')
fi;

if divs.numer(divs.divs(-1, 3)) <> -1 then
   print('Error in divs.numer, 1')
fi;

if divs.numer(divs.divs(-1:3)) <> -1 then
   print('Error in divs.numer, 2')
fi;

if divs.divs(2, 3) ** (-2) <> divs.divs(9, 4) then
   print('Error in divs/**')
fi;

if int divs.divs(-3, -2) <> 1 then
   print('Error in divs/int')
fi;


print("\nChecking bit and logical operators\n");

to 10 do
   if 4 >>> 1 <> 2 then
      print('Error in >>>')
   fi;
   if 4 <<< 1 <> 8 then
      print('Error in <<<, 1')
   fi;
   if 4 shift 1 <> 8 then
      print('Error in shift, 1')
   fi;
   if 4 shift -1 <> 2 then
      print('Error in shift, 2')
   fi;
   if 2 && 0 <> 0 then
      print('Error in &&')
   fi;
   if 2 || 1 <> 3 then
      print('Error in ||')
   fi;
   if 2 xor 1 <> 1 then
      print('Error in xor')
   fi;
   if 1 <<< 2 >>> 1 <> 2 then
      print('Error in <<< and >>>')
   fi;
   if (1 <<< 31) <> 2^31 then
      print('Error in <<<, 2')
   fi;
   if not(1 <<< 32 = 2^0 or 1 <<< 32 = 0) then  # 0 for PowerPC, ARM
      print('Error in <<<, 3')
   fi
od;

if 5703180337 <<< 10 <> 3196109824 then  # overflow situation, PPC: 4294967295
   print('Error in <<<, 4', 5703180337 <<< 10, 'instead of 3196109824', 'signedbits=' & tostring(environ.kernel().signedbits))
fi;


# advanced Boolean operators
operands := [[false, false, true, true], [false, true, false, true]];

r := ['nor'  ~ [true, false, false, false],
      'nand' ~ [true, true, true, false],
      'xor'  ~ [false, true, true, false],
      'xnor' ~ [true, false, false, true]];

f := ['nor'  ~  << x, y -> x nor y >>,
      'nand' ~ << x, y -> x nand y >>,
      'xor'  ~  << x, y -> x xor y >>,
      'xnor' ~ << x, y -> x xnor y >>];

for i in tables.indices(r) do
   for j to 4 do
      if f[i](operands[1, j], operands[2, j]) <> r[i, j] then
         print('Error in ' & i & '/' & j)
      fi
   od
od;


# bitwise operations

import bytes;

environ.kernel(signedbits = false);

hashmask := << n -> (1 <<< n) - 1 >>

doprint := false;

oaat := proc(s :: string, n :: number) is
   local h := 0;
   for i in s do
      inc h, abs i;
      h := bytes.numto32(h);  # equals: h > 2^32 - 1 ? h := h >>> 31;
      doprint ? print('loop ' & i & ': h = ' & h);
      inc h, h <<< 10;
      doprint ? print(' 1  h = ' & h);
      h := bytes.numto32(h);
      doprint ? print(' 2  h = ' & h);
      h := h ^^ (h >>> 6)
      doprint ? print(' 3  h = ' & h);
      h := bytes.numto32(h);
      doprint ? print(' 4  h = ' & h);
   od;
   inc h, h <<< 3;
   doprint ? print('5 h = ' & h);
   h := bytes.numto32(h);
   doprint ? print('6 h = ' & h);
   h := h ^^ (h >>> 11);
   doprint ? print('7 h = ' & h);
   h := bytes.numto32(h);
   doprint ? print('8 h = ' & h);
   inc h, h <<< 15;
   doprint ? print('9 h = ' & h);
   h := bytes.numto32(h);
   doprint ? print('A h = ' & h);
   return if n <> 0 then h && hashmask(n) else h fi
end;

/*
loop a: h = 97
 1  h = 99425
 2  h = 99425
 3  h = 98928
 4  h = 98928
loop l: h = 99036
 1  h = 101511900
 2  h = 101511900
 3  h = 102023447
 4  h = 102023447
loop e: h = 102023548
 1  h = 1494921596
 2  h = 1494921596
 3  h = 1484708793
 4  h = 1484708793
loop x: h = 1484708913
 1  h = 5703180337
 2  h = 1408213041
 3  h = 1386225329
 4  h = 1386225329
loop a: h = 1386225426
 1  h = 3541853970
 2  h = 3541853970
 3  h = 3494909582
 4  h = 3494909582
loop n: h = 3494909692
 1  h = 4574676732
 2  h = 279709436
 3  h = 284079847
 4  h = 284079847
loop d: h = 284079947
 1  h = 3419136843
 2  h = 3419136843
 3  h = 3370437830
 4  h = 3370437830
loop e: h = 3370437931
 1  h = 5840140587
 2  h = 1545173291
 3  h = 1567169743
 4  h = 1567169743
loop r: h = 1567169857
 1  h = 4326302017
 2  h = 31334721
 3  h = 31021508
 4  h = 31021508
loop w: h = 31021627
 1  h = 1732396603
 2  h = 1732396603
 3  h = 1725910819
 4  h = 1725910819
loop a: h = 1725910916
 1  h = 3827130244
 2  h = 3827130244
 3  h = 3884788474
 4  h = 3884788474
loop l: h = 3884788582
 1  h = 4768580454
 2  h = 473613158
 3  h = 474622059
 4  h = 474622059
loop z: h = 474622181
 1  h = 1156431077
 2  h = 1156431077
 3  h = 1174276630
 4  h = 1174276630
5 h = 1978555078
6 h = 1978555078
7 h = 1977673485
8 h = 1977673485
9 h = 3915867917
A h = 3915867917
3915867917 */

#   lua_pushnumber(L, (n != 0) ? h & hashmask(n) : h);  /* 2.11.4 change, 2.12.0 RC 4 patch: use (n != 0) instead of just (n) */
if oaat('alexanderwalz', 0) <> 3915867917 then
   print('Error in left/right-shifts/bitwise (X)OR: ' & oaat('alexanderwalz', 0) & ' instead of 3915867917')
fi;

environ.kernel(signedbits = false);
if (1 <<< 31) <> 2^31 then
   print('Error in left-shift: ' & 1 <<< 31 & ' ' & 2^31)
fi;

if (1024 >>> 9) <> 2 then
   print('Error in right-shift')
fi;

environ.kernel(signedbits = true);

# math.koadd
x, q -> 0;
y := 0.1;
while x < 1 do
   x, q := math.koadd(x, y, q)
od;
if x <> 1 then
   print('Error in math.koadd')
fi;

# Misc, we have already loaded the `bytes' package

import fastmath

if bytes.tobytes(Pi) <> seq(24, 45, 68, 84, 251, 33, 9, 64) then
   print('Error in bytes.tobytes')
fi;

if bytes.tonumber(seq(24, 45, 68, 84, 251, 33, 9, 64)) <> Pi then
   print('Error in bytes.tonumber')
fi;

if math.dms(2.9775) <> 2.5839 then
   print('Error in math.dms')
fi

if math.dd(2.5839) <> 2.9775 then
   print('Error in math.dd')
fi;

zerosignbit := math.signbit(0);

if zerosignbit <> false then
   print('Error in math.signbit/1')
fi;

minuszerosignbit := math.signbit(-0);

if minuszerosignbit <> true then
   print('Error in math.signbit/2')
fi;

pluszerosignbit := math.signbit(+0);

if pluszerosignbit <> false then
   print('Error in math.signbit/3')
fi;

if math.convertbase('FF', 16, 10) <> '255' or math.convertbase('FF', 16, 2) <> '11111111' then
   print('Error in math.convertbase/1')
fi;

# check stack validity

stack.resetd();

pushd(100, 101);

math.convertbase('FF', 16, 2);

bxxx := stack.dumpd()

if bxxx <> seq(100, 101) then
   print('Error in math.convertbase/2', bxxx)
fi;

clear bxxx;

stack.resetd()

if math.flipsign(-1, -1) <> 1 then
   print('Error in math.flipsign/1')
fi;
if math.flipsign(-1, 1) <> -1 then
   print('Error in math.flipsign/2')
fi;
if math.flipsign(1, 1) <> 1 then
   print('Error in math.flipsign/3')
fi;
if math.flipsign(1, -1) <> -1 then
   print('Error in math.flipsign/4')
fi;
if math.fld(3, 2) <> 1 then
   print('Error in math.fld/1')
fi;
if math.fld(-3, 2) <> -2 then
   print('Error in math.fld/2')
fi;
if math.cld(3, 2) <> 2 then
   print('Error in math.cld/1')
fi;
if math.cld(-3, 2) <> -1 then
   print('Error in math.cld/2')
fi;
if math.isqrt(3) <> 1 then
   print('Error in math.isqrt/1')
fi;
if math.isqrt(4) <> 2 then
   print('Error in math.isqrt/2')
fi;
if math.isqrt(5) <> 2 then
   print('Error in math.isqrt/3')
fi;
for n from 0 to 10 do
   if not approx(lngamma(n + 1), math.lnfact(n)) then
      print('Error in math.lnfact/' & n)
   fi
od;

if infinite(1) <> false then
   print('Error in infinite/1')
fi;
if infinite(undefined) <> false then
   print('Error in infinite/2')
fi;
if infinite(infinity) <> true then
   print('Error in infinite/3')
fi;
if infinite(sqrt(-1)) <> false then
   print('Error in infinite/4')
fi;
if infinite(-infinity) <> true then
   print('Error in infinite/5')
fi;

for i from -10 to 10 by 0.1 do
   if [frexp(i*2^i)] <> [math.mantissa(i*2^i), math.exponent(i*2^i)] then
      print('Error in math.mantissa/math.exponent', i)
   fi
od;

# base 0, 1, 3 and 4 have been tested, as well
for i from 0 to 20 do
   if not math.ispow2(2^i) then
      print('Error in math.ispow2', i)
   fi
od;

for i from -10 to 10 by 0.1 do
   if math.eps(i, true) <> math.epsilon(i) then
      # print('Error in math.epsilon', i, math.eps(i, true) - math.epsilon(i))
   fi
od;

if not math.issubnormal(math.nextafter(0, 1)) then
   print('Error in math.issubnormal, 1')
fi;

if not math.issubnormal(math.nextafter(2^-1022, 0)) then
   print('Error in math.issubnormal, 2')
fi;

if math.zerosubnormal(math.nextafter(2^-1022, 0)) <> 0 then
   print('Error in math.zerosubnormal')
fi;

# Maple V R4: with(combinat, fibonacci); seq(fibonacci(i), i=0 .. 76);
ffib := [ 0 ~ 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181,
       6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269,
       2178309, 3524578, 5702887, 9227465, 14930352, 24157817, 39088169, 63245986, 102334155,
       165580141, 267914296, 433494437, 701408733, 1134903170, 1836311903, 2971215073, 4807526976,
       7778742049, 12586269025, 20365011074, 32951280099, 53316291173, 86267571272, 139583862445,
       225851433717, 365435296162, 591286729879, 956722026041, 1548008755920, 2504730781961,
       4052739537881, 6557470319842, 10610209857723, 17167680177565, 27777890035288, 44945570212853,
       72723460248141, 117669030460994, 190392490709135, 308061521170129, 498454011879264,
       806515533049393, 1304969544928657, 2111485077978050, 3416454622906707 ];

for i from 0 to 76 do
   if ffib[i] <> math.fib(i) then print('Error in math.fib, ', i, ffib[i], math.fib(i), ffib[i] - math.fib(i)) fi
od;

clear ffib;

if isnonposint(1.1) <> false then
   print('error in isnonposint')
fi;

if sign(1) <> 1 then
   print('error in sign, 1')
fi;
if sign(-1) <> -1 then
   print('error in sign, 2')
fi;
if sign(0) <> 0 then
   print('error in sign, 3')
fi;
if sign(undefined) <> undefined then
   print('error in sign, 4')
fi;

if ilog2(1) <> 0 then
   print('error in ilog2, 1')
fi;
if ilog2(1G) <> 30 then
   print('error in ilog2, 2')
fi;


for i to 50 by 0.1 do
   if approx(i/%51000, i/51000*100) = false then
      print('Error in /% operator ' & i)
   fi
od;

for i to 50 by 0.1 do
   if approx(i*%51000, i*51000/100) = false then
      print('Error in *% operator ' & i)
   fi
od;

to 10 do
   if 100 *% 2 /% 3 +% 4 -% 5 ~<> 65.866666666667 then
      print('Error in percentage or ~<> operators')
   fi
od;

for i from -10 to 10 by 0.1 do
   if not approx(cbrt(i), root(i, 3)) then
      print('Error in cbrt with ' & i)
   fi
od;

for i from -10 to 10 by 0.1 do
   for j from -10 to 10 do
      if not approx((i!j)^(1/3), cbrt(i!j)) then
         print('Error in cbrt with complex ', i!j)
      fi
   od
od;

for i from -10 to 10 by 0.1 do
   for n to 10 do
      if i < 0 and not(even(n)) then skip fi;
      if not approx((i)^(1/n), proot(i, n)) then
         print('Error in proot with real ', i, n)
      fi
   od
od;

for i from -10 to 10 by 0.1 do
   for n to 10 do
      if not approx((i!0)^(1/n), proot(i+I*0, n)) then
         print('Error in proot with complex ', i!0, n, (i!0)^(1/n), proot(i+I*0, n))
      fi
   od
od;

if approx(drem(6.5, 2.3), -0.4) <> true then
   print('Error in drem')
fi;

if math.norm(5, 0:10, -100:100) <> 0 then  /* 2.10.1 fix */
   print('Error in math.norm')
fi;

z := 1!-2;

if flip(z) <> imag(z)!real(z) then
   print('Error in flip')
fi;

if conjugate(cos(z)) <> cosxx(z) then
   print('Error in cosxx')
fi;

if cabs(z) <> abs(real(z))!abs(imag(z)) then
   print('Error in cabs')
fi;

for i from -10 to 10 by 0.1 do
   for j from -10 to 10 by 0.1 do
      if abs((i+I*j) - (j+I*i)) <> ((i+I*j) |- (j+I*i)) then
         print('Error in |-', i, j, abs((i+I*j) - (j+I*i)), (i+I*j) |- (j+I*i))
      fi
   od
od;

# absiff
for i from -50 to 50 do
   a := i |- (i - 1);
   b := abs(i - (i - 1));
   if a <> b then
      print('Error in |-, ' & i)
   fi
od;

r := 2.8284271247462;
if abs(1!2 - 3!4) ~<> r or 1!2 |- 3!4 ~<> r or 1+2*I |- 3!4 ~<> r or 1+2*I |- 3+4*I ~<> r then
   print('Error in abs/|- precedence')
fi;

# numeric compare
for i from -50 to 50 do
   if i < 0 and i | 0 <> -1 then
      print('Error in |, negative, ' & i)
   elif i = 0 and i | 0 <> 0 then
      print('Error in |, zero, ' & i)
   elif i > 0 and i | 0 <> 1 then
      print('Error in |, positive, ' & i)
   fi
od;

if not(math.fpclassify(0) = math.fp_zero) then
   print('Error in math.fpclassify, 1')
fi;

if not(math.fpclassify(undefined) = math.fp_nan) then
   print('Error in math.fpclassify, 2')
fi;

if not(math.fpclassify(1) = math.fp_normal) then
   print('Error in math.fpclassify, 3')
fi;

if not(math.fpclassify(infinity) = math.fp_infinite) then
   print('Error in math.fpclassify, 4')
fi;

if not(math.fpclassify(-infinity) = math.fp_infinite) then
   print('Error in math.fpclassify, 5')
fi;

if math.dirac(0) <> infinity or math.dirac(++0) <> 0 or math.dirac(0, Eps) = infinity then
   print('Error in math.dirac')
fi;

if math.clip(1.1) <> 1 or math.clip(1) <> 1 or math.clip(0) <> 0 or math.clip(-1.1) <> -1 then
   print('Error in math.clip')
fi;

if math.unitise(++0) <> 0 or math.unitise(0) <> 0 or math.unitise(Eps) <> 0 or math.unitise(++Eps) <> 1 then
   print('Error in math.unitise')
fi;

if math.rectangular(-0.5) <> 0.5 or math.rectangular(-0.6) <> 0 or math.rectangular(0.5) <> 0.5 or
   math.rectangular(0) <> 1 or math.rectangular(0.4) <> 1 or math.rectangular(0.6) <> 0 then
      print('Error in math.rectangular, 1')
fi;

if math.rectangular(0, true) <> 1 or math.rectangular(0.5, true) <> undefined or math.rectangular(0.25, true) <> 1 or
   math.rectangular(0.75, true) <> 0 or math.rectangular(-0.5, true) <> undefined or math.rectangular(-0.25, true) <> 1 or
   math.rectangular(-0.75, true) <> 0 then
      print('Error in math.rectangular, 2')
fi;

if math.triangular(0.5) <> 0 or math.triangular(0.25) <> 0.5 or math.triangular(0) <> 1 or
   math.triangular(-0.25) <> 0.5 or math.triangular(-0.5) <> 0 or math.triangular(-1) <> 0 or
   math.triangular(1) <> 0 or math.triangular(1, 0, 2) <> 1 or math.triangular(0.5, 0, 2) <> 0.5 then
      print('Error in math.triangular')
fi;

for i from -2 to 2 by 0.1 do
   if i = 0 then
      if math.uexponent(i) <> -1023 then
         print('Error in math.uexponent at x=', i)
      fi
   else
      if math.uexponent(i) <> math.exponent(i) - 1 then
         print('Error in math.uexponent at x=', i)
      fi
   fi
od

a, b := bytes.leadzeros(15);
if a <> 28 or b <> 15 then
   print('Error in math.leadzeros')
fi;

if bytes.onebits(15) <> 4 then
   print('Error in math.onebits')
fi;

a := bytes.cast(Pi);

if bytes.getlow(a) <> 1413754136 or bytes.gethigh(a) <> 1074340347 then
   print('Error in bytes.gethigh or getlow')
fi;

a := bytes.cast(0);

if bytes.sethigh(a, 1074340347) ~<> 3.1415920257568 or bytes.setlow(a, 1413754136) ~<> Pi then
   print('Error in bytes.sethigh or setlow')
fi;

if bytes.setwords(a, 1074340347, 1413754136) ~<> Pi then
   print('Error in bytes.setwords')
fi;

x0, y0 := bytes.setdouble(a, Pi);
if x0 <> bytes.gethigh(a) or y0 <> bytes.getlow(a) then
   print('Error in bytes.setdouble')
fi;

if bytes.getdouble(a) ~<> Pi then
   print('Error in bytes.getdouble')
fi;

x1, y1 := bytes.getwords(a);

if x1 <> x0 or y1 <> y0 or x1 <> bytes.gethigh(a) or y1 <> bytes.getlow(a) then
   print('Error in bytes.getwords')
fi;

if bytes.setnumwords(1074340347, 1413754136) ~<> Pi then
   print('Error in bytes.setnumwords')
fi;

x0, y0 := bytes.numwords(Pi);
if x0 <> 1074340347 or y0 <> 1413754136 then
   print('Error in bytes.numwords')
fi;

x := bytes.setnumhigh(0, bytes.numhigh(Pi));  # = 1074340347
if x ~<> 3.1415920257568 or bytes.setnumlow(x, bytes.numlow(Pi)) ~<> Pi then  # = 1413754136
   print('Error in bytes.setnumhigh/low')
fi;

if bytes.setnumlow(bytes.setnumhigh(0, 3221823995), 1413754136) <> - Pi then
   print('error in `bytes.setnumhigh/low`, 1')
fi;

if bytes.setwords(bytes.cast(), 3221823995, 1413754136) <> -Pi then
   print('error in `bytes.setwords`, 1')
fi;

x := bytes.sethigh(bytes.cast(0), 3221823995);
if bytes.setlow(bytes.cast(x), 1413754136) <> - Pi then
   print('error in `bytes.sethigh/low`, 1')
fi;

x := bytes.cast(0);
if bytes.getunbiased(x) <> -1023 then
   print('error in `bytes.getunbiased`')
fi;

# For more examples, see: github.com/Shopify/lua-tests/blob/master/bitwise.lua

environ.kernel(signedbits = false);

if bytes.shift32(1, -31) <> 2^31 then
   print('Error in bytes.shift32, 1')
fi;

if bytes.shift32(1, 0) <> 1 then
   print('Error in bytes.shift32, 2')
fi;

if bytes.shift32(2, 1) <> 1 then
   print('Error in bytes.shift32, 3')
fi;

if bytes.shift32(0x12345678, -4) <> 0x23456780 then
   print('Error in bytes.shift32, 4')
fi;

if bytes.shift32(0x12345678, 4) <> 0x01234567 then
   print('Error in bytes.shift32, 5')
fi;

if bytes.arshift32(0x12345678, 0) <> 0x12345678 then
   print('Error in bytes.arshift32, 1')
fi;

if bytes.arshift32(0x12345678, 1) <> 0x12345678/2 then
   print('Error in bytes.arshift32, 2')
fi;

if bytes.arshift32(0x12345678, -1) <> 0x12345678*2 then
   print('Error in bytes.arshift32, 3')
fi;

if bytes.rotate32(1, -31) <> 2^31 then
   print('Error in bytes.rotate32, 1')
fi;

if bytes.rotate32(1, 1) <> 2^31 then
   print('Error in bytes.rotate32, 2')
fi;

if bytes.rotate32(1, 0) <> 1 then
   print('Error in bytes.rotate32, 3')
fi;

if bytes.rotate32(2, 1) <> 1 then
   print('Error in bytes.rotate32, 4')
fi;

if bytes.extract32(0b011101, 2, 3) <> 0b111 then
   print('Error in bytes.extract32, 1')
fi;

if bytes.extract32(0b011101, 0, 3) <> 0b101 then
   print('Error in bytes.extract32, 2')
fi;

if bytes.replace32(0b011101, 0b1, 1, 1) <> 0b11111 then
   print('Error in bytes.replace32, 1')
fi;

if bytes.arshift32(0b10111, -1) <> 0b101110 then
   print('Error in bytes.arshift32, 1')
fi;

if bytes.arshift32(0b10111, 1) <> 0b1011 then
   print('Error in bytes.arshift32, 2')
fi;


fastmath.sincosfast(Pi);
fastmath.sqroot(2);
fastmath.sqrtfast(2);

for i from -10 to 10 by 0.1 do
   if square(i) <> i^2 then
      print('Error in `square with`', i)
   fi;
   if square(i!(i+1)) ~<> (i!(i+1))^2 then
      print('Error in `square` with complex', i, square(i!(i+1)), (i!(i+1))^2)
   fi
od;

for i in io.lines('mathmul.tst') do
   local a, f, nofargs, nofrets, argseq, maple;
   skip when i and i[1] = '#';
   break when unassigned i;
   a := tonumber @ (i split ';');
   f := a[1];        # function name
   nofargs := a[2];  # number of arguments
   nofrets := a[3];  # number of results
   argseq := '';
   for j to nofargs do
      argseq := argseq & a[3 + j] & ','
   od;
   argseq := strings.chomp(argseq, ',');
   maple := [];
   for j to nofrets do
      insert tonumber(a[3 + nofargs + j]) into maple
   od;
   s := f & '(' & argseq & ')';  # build Agena call
   loadstring('agena := [' & s & ']')();
   if false and agena = [undefined] and a[-1] = 'no' then
      s := prefix & a[1] & '(' & a[2] & '+0*I, ' & a[3] & '+0*I)';  # build Agena call
      loadstring('agena := ' & s)();
   fi;
   if maple[1] in {'singularityencountered', 'divisionbyzero'} then maple[1] := undefined fi;
   for j to nofrets do
      if isdiff(agena[j], maple[j]) then
         # with very large values, roundoff errors are commonplace, so do some raw comparisons
         if tostring(agena[j]) <> tostring(maple[j]) and (not(agena[j]=undefined and type(maple[j]) = complex)) then
            print('>>> ', s, agena[j], maple[j])
         fi
      fi
   od;
od;


for i from -1k to 1k by 0.01 do
   z := i!(i-1);
   if z**2 - z ~<> z squareadd -z then print('Error in squareadd,', i, (z**2 - z) |- (z squareadd -z)) fi
od;

for i from -1k to 1k by 0.01 do
   z := i;
   if z**2 - z ~<> z squareadd -z then print('Error in squareadd,', i, (z**2 - z) |- (z squareadd -z)) fi
od;


for i from -1k to 1k by 0.01 do
   z := i!(i-1);
   if 1/sqrt z ~<> invsqrt z then print('Error in invsqrt,', i, 1/sqrt z |- invsqrt z) fi
od;

for i from -1k to 1k by 0.01 do
   z := i;
   if 1/sqrt z ~<> invsqrt z then print('Error in invsqrt,', i, 1/sqrt z |- invsqrt z) fi
od;


for i from -1k to 1k by 0.01 do
   x := i!(i-1);
   y := (i+1)!i;
   if abs(x - y) <> x |- y then print('Error in |-,', i, abs(x - y), x |- y) fi
od;

for i from -1k to 1k by 0.01 do
   x := i + 1;
   y := i - 1;
   if abs(x - y) <> x |- y then print('Error in |-,', i, abs(x - y), x |- y) fi
od;



print('math: ' & tostring(time() - _time) & ' seconds');

# delete all global values
clear _time, a, A, A1, A2, A3, A4, A5, A6, A7, M, stattest, agena, b, c,
   calcfns, statfns, cmplxtest, cmplxtest2, f, g, i, isdiff, lb, lew, m, n,
   p, q, r, s1, statfns, j, maple, mathall, mathfns, mathops, operands,
   prefix, s, s1, se, v, x, x0, x1, y, y0, y1, z, xx, xy, xz, hashmask, oaat, zerosignbit,
   pluszerosignbit, minuszerosignbit;

