octave_packages/fixed-0.7.10/fixedpoint.m

   1 ## Copyright (C) 2003  Motorola Inc
   2 ## Copyright (C) 2003  David Bateman
   3 ##
   4 ## This program is free software; you can redistribute it and/or modify
   5 ## it under the terms of the GNU General Public License as published by
   6 ## the Free Software Foundation; either version 2 of the License, or
   7 ## (at your option) any later version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful,
  10 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 ## GNU General Public License for more details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License
  15 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  16
  17 ## -*- texinfo -*-
  18 ## @deftypefn {Function File} {} fixedpoint ('help')
  19 ## @deftypefnx {Function File} {} fixedpoint ('info')
  20 ## @deftypefnx {Function File} {} fixedpoint ('info', @var{mod})
  21 ## @deftypefnx {Function File} {} fixedpoint ('test')
  22 ## @deftypefnx {Function File} {} fixedpoint ('test', @var{mod})
  23 ##
  24 ## Manual and test code for the Octave Fixed Point toolbox. There are
  25 ## 5 possible ways to call this function.
  26 ##
  27 ## @table @code
  28 ## @item fixedpoint ('help')
  29 ## Display this help message. Called with no arguments, this function also
  30 ## displays this help message
  31 ## @item fixedpoint ('info')
  32 ## Open the Fixed Point toolbox manual
  33 ## @item fixedpoint ('info', @var{mod})
  34 ## Open the Fixed Point toolbox manual at the section specified by
  35 ## @var{mod}
  36 ## @item fixedpoint ('test')
  37 ## Run all of the test code for the Fixed Point toolbox.
  38 ## @var{mod}.
  39 ## @end table
  40 ##
  41 ## Valid values for the varibale @var{mod} are
  42 ##
  43 ## @table @asis
  44 ## @item 'basics'
  45 ## The section describing the use of the fixed point toolbox within Octave
  46 ## @item 'programming'
  47 ## The section descrining how to use the fixed-point type with oct-files
  48 ## @item 'example'
  49 ## The section describing an in-depth example of the use of the fixed-point
  50 ## type
  51 ## @item 'reference'
  52 ## The refernce section of all of the specific fixed point operators and
  53 ## functions
  54 ## @end table
  55 ##
  56 ## Please note that this function file should be used as an example of the
  57 ## use of this toolbox.
  58 ## @end deftypefn
  59
  60 function retval = fixedpoint(typ, tests)
  61
  62   if (nargin < 1)
  63     typ = "help";
  64     tests = "all";
  65   elseif (nargin < 2)
  66     tests = "all";
  67   endif
  68
  69   infofile = [fileparts(mfilename('fullpath')),'/doc.info'];
  70   n = 10;
  71   m = 20;
  72   is = 7;
  73   ds = 6;
  74   vals = [-10,10];
  75
  76   vers = sscanf(OCTAVE_VERSION, "%i.%i.%i");
  77
  78   if strcmp(tests,"all")
  79     nodename = "Top";
  80   elseif strcmp(tests,"basics")
  81     nodename = "Basics";
  82   elseif  strcmp(tests,"programming")
  83     nodename = "Programming";
  84   elseif  strcmp(tests,"example")
  85     nodename = "Example";
  86   elseif  strcmp(tests,"reference")
  87     nodename = "Function Reference";
  88   else
  89     error ("fixedpoint: unrecognized section");
  90   endif
  91
  92   if (strcmp(typ,"help"))
  93     help ("fixedpoint");
  94   elseif (strcmp(typ,"info"))
  95     infopaths = ["."];
  96     if (!isempty(char (strsplit (path, ":"))))
  97       infopaths =[infopaths; char (strsplit (path, ":"))];
  98     endif
  99     if (!isempty(char (strsplit (DEFAULT_LOADPATH, ":"))))
 100       infopaths =[infopaths; char (strsplit (DEFAULT_LOADPATH, ":"))];
 101     endif
 102     for i=1:size(infopaths,1)
 103       infopath = deblank(infopaths(i,:));
 104       len = length(infopath);
 105       if (len)
 106         if (len > 1 && strcmp(infopath([len-1, len]),"//"))
 107           [status, showfile] = system(["find '", infopath(1:len-1), ...
 108                                        "' -name ", infofile]);
 109         else
 110           [status, showfile] = system(["find '", infopath, "' -name ", ...
 111                                        infofile, " -maxdepth 1"]);
 112         endif
 113         if (length(showfile))
 114           break;
 115         endif
 116       endif
 117     end
 118     if (!exist("showfile") || !length(showfile))
 119       error("fixedpoint: info file not found");
 120     endif
 121     if (showfile(length(showfile)) == "\n")
 122       showfile = showfile(1:length(showfile)-1);
 123     endif
 124
 125     if (exist("INFO_PROGAM"))
 126       [testret, testout] = system(["'", INFO_PROGRAM, "' --version"]);
 127       if (testret)
 128         error("fixedpoint: info command not found");
 129       else
 130         system(["'", INFO_PROGRAM, "' --file '", showfile, "' --node '", ...
 131                 nodename, "'"]);
 132       endif
 133     else
 134       [testret, testout] = system("info --version");
 135       if (testret)
 136         error("fixedpoint: info command not found");
 137       else
 138         system(["info --file '", showfile, "' --node '", nodename, "'"]);
 139       endif
 140     endif
 141   elseif (strcmp(typ,"test"))
 142     pso = page_screen_output();
 143     unwind_protect
 144       page_screen_output(0);
 145       feps = 1 / 2 ^ ds;
 146
 147       fprintf("\n<< Fixed Point Load Type >>\n");
 148
 149       ## Load a fixed point variable to start things off
 150       x = fixed(is,ds);
 151
 152       if (!exist("fixed_point_warn_overflow") || !exist("fixed_point_debug") ||
 153           !exist("fixed_point_count_operations") ||
 154           !exist("fixed_point_version"))
 155         error("Can not find fixed point type")
 156       endif
 157
 158       fprintf("  Found Fixed Point Toolbox (version %s)\n",
 159               fixed_point_version);
 160
 161       fprintf("\n<< Fixed Point Creation >>\n");
 162       fprintf("  Scalar Creation:                          ");
 163       zero = fixed(is,ds);
 164       if (zero.x != 0.) || (zero.int != is) || (zero.dec != ds)
 165         error ("FAILED");
 166       endif
 167       if (!isscalar(zero) || !ismatrix(zero) || !isfixed(zero) ||
 168           !isvector(zero))
 169         error ("FAILED");
 170       endif
 171       if any(zero) || all(zero)
 172         error ("FAILED");
 173       endif
 174       if (size(zero) != [1,1]) || (length(zero) != 1)
 175         error ("FAILED");
 176       endif
 177       zero.int = zero.int+1;
 178       zero.dec = zero.dec+1;
 179       if (zero.x != 0.) || (zero.int != is+1) || (zero.dec != ds+1)
 180         error ("FAILED");
 181       endif
 182       empty = fixed(is,ds,[]);
 183       if (isscalar(empty) || ismatrix(empty) || !isfixed(empty) ||
 184           isvector(empty) || !isempty(empty) || isempty(zero))
 185         error ("FAILED");
 186       endif
 187       for i=1:100
 188         x = (vals(2)-vals(1))*rand(1,1)+vals(1);
 189         a = fixed(is, ds, x);
 190         if (abs(a.x - x) > feps) || (a.sign != sign(a.x)) || !isfixed(a)
 191           error ("FAILED");
 192         endif
 193       endfor
 194       fprintf("PASSED\n");
 195       fprintf("  Matrix Creation:                          ");
 196       zero = fixed(is,ds,zeros(n,n));
 197
 198       if (any(any(zero.x)) || any(any(zero.int != is)) ||
 199           any(any(zero.dec != ds)))
 200         error ("FAILED");
 201       endif
 202       if (isscalar(zero) || !ismatrix(zero) || !isfixed(zero) ||
 203           isvector(zero))
 204         error ("FAILED");
 205       endif
 206       if any(any(zero)) || all(all(zero))
 207         error ("FAILED");
 208       endif
 209       if (size(zero) != [n,n])
 210         error ("FAILED");
 211       endif
 212       zero.int = zero.int+1;
 213       zero.dec = zero.dec+1;
 214       if (any(any(zero.x)) || any(any(zero.int != is+1)) ||
 215           any(any(zero.dec != ds+1)))
 216         error ("FAILED");
 217       endif
 218       for i=1:100
 219         x = (vals(2)-vals(1))*rand(n,n)+vals(1);
 220         a = fixed(is, ds, x);
 221         if (any(any(abs(a.x - x) > feps)) || any(any(a.sign != sign(a.x)))
 222             || !isfixed(a))
 223           error ("FAILED");
 224         endif
 225       endfor
 226       b = freshape(a,n*n,1);
 227       if (isscalar(b) || !ismatrix(b) || !isfixed(b) || !isvector(b))
 228         error ("FAILED");
 229       endif
 230       if any(b != a(:))
 231         error ("FAILED");
 232       endif
 233       zero = fixed(is,ds,0);
 234       vec = fixed(is,ds,[1:n]);
 235       a = fdiag(vec);
 236       if (size(a,1) != n || size(a,2) != n)
 237         error("FAILED");
 238       endif
 239       for i=1:n
 240         for j=1:n
 241           if ((i == j) && (a(i,j) != vec(i)))
 242             error("FAILED");
 243           elseif ((i != j) && (a(i,j) != zero))
 244             error("FAILED");
 245           endif
 246         end
 247       end
 248       vec = fdiag(a);
 249       if (length(vec) != n)
 250         error("FAILED");
 251       endif
 252       for i=1:n
 253         if (a(i,i) != vec(i))
 254           error("FAILED");
 255         endif
 256       end
 257       fprintf("PASSED\n");
 258       fprintf("  Complex Scalar Creation:                  ");
 259       onei = fixed(is*(1+1i),ds*(1+1i),1+1i);
 260       if ((onei.x != 1+1i) || (onei.int != is*(1+1i)) ||
 261           (onei.dec != ds*(1+1i)))
 262         error ("FAILED");
 263       endif
 264       if (!isscalar(onei) || !ismatrix(onei) || !isfixed(onei) ||
 265           !isvector(onei) || !iscomplex(onei))
 266         error ("FAILED");
 267       endif
 268       if !any(onei) || !all(onei)
 269         error ("FAILED");
 270       endif
 271       if (size(onei) != [1,1]) || (length(onei) != 1)
 272         error ("FAILED");
 273       endif
 274       onei.int = onei.int+1;
 275       onei.dec = onei.dec+1;
 276       if ((onei.x != 1+1i) || (onei.int != is*(1+1i)+1) ||
 277           (onei.dec != ds*(1+1i)+1))
 278         error ("FAILED");
 279       endif
 280       empty = fixed(is*(1+1i),ds*(1+1i),[]);
 281       if (isscalar(empty) || ismatrix(empty) || !isfixed(empty) ||
 282           isvector(empty) || !isempty(empty) || isempty(onei) ||
 283           iscomplex(empty))
 284         ## Not complex, since its type is narrowed!!
 285         error ("FAILED");
 286       endif
 287
 288       for i=1:100
 289         x = (vals(2)-vals(1))*rand(1,1)+vals(1) + ...
 290             1i * (vals(2)-vals(1))*rand(1,1)+vals(1);
 291         a = fixed(is, ds, x);
 292         if ((abs(real(a.x) - real(x)) > feps) ||
 293             (abs(imag(a.x) - imag(x)) > feps) || !isfixed(a))
 294           error ("FAILED");
 295         endif
 296       endfor
 297       fprintf("PASSED\n");
 298       fprintf("  Complex Matrix Creation:                  ");
 299       onei = fixed(is*(1+1i),ds*(1+1i),ones(n,n)*(1+1i));
 300
 301       if (any(any(onei.x != 1+1i)) || any(any(onei.int != is*(1+1i))) ||
 302           any(any(onei.dec != ds*(1+1i))))
 303         error ("FAILED");
 304       endif
 305       if (isscalar(onei) || !ismatrix(onei) || !isfixed(onei) ||
 306           isvector(onei) || !iscomplex(onei))
 307         error ("FAILED");
 308       endif
 309       if !any(any(onei)) || !all(all(onei))
 310         error ("FAILED");
 311       endif
 312       if (size(onei) != [n,n])
 313         error ("FAILED");
 314       endif
 315       onei.int = onei.int+1;
 316       onei.dec = onei.dec+1;
 317       if (any(any(onei.x != 1+1i)) || any(any(onei.int != is*(1+1i)+1)) ||
 318           any(any(onei.dec != ds*(1+1i)+1)))
 319         error ("FAILED");
 320       endif
 321       for i=1:100
 322         x = (vals(2)-vals(1))*rand(n,n)+vals(1) + ...
 323             1i * (vals(2)-vals(1))*rand(n,n)+vals(1);
 324         a = fixed(is, ds, x);
 325         if (any(any(abs(real(a.x) - real(x)) > feps)) ||
 326             any(any(abs(imag(a.x) - imag(x)) > feps)) || !isfixed(a))
 327           error ("FAILED");
 328         endif
 329       endfor
 330       b = freshape(a,n*n,1);
 331       if (isscalar(b) || !ismatrix(b) || !isfixed(b) || !isvector(b))
 332         error ("FAILED");
 333       endif
 334       if any(b != a(:))
 335         error ("FAILED");
 336       endif
 337       zero = fixed(is,ds,0);
 338       vec = fixed(is*(1+1i),ds*(1+1i),[1:n]*(1+1i));
 339       a = fdiag(vec);
 340       if (size(a,1) != n || size(a,2) != n)
 341         error("FAILED");
 342       endif
 343       for i=1:n
 344         for j=1:n
 345           if ((i == j) && (a(i,j) != vec(i)))
 346             error("FAILED");
 347           elseif ((i != j) && (a(i,j) != zero))
 348             error("FAILED");
 349           endif
 350         end
 351       end
 352       vec = fdiag(a);
 353       if (length(vec) != n)
 354         error("FAILED");
 355       endif
 356       for i=1:n
 357         if (a(i,i) != vec(i))
 358           error("FAILED");
 359         endif
 360       end
 361       fprintf("PASSED\n");
 362       fprintf("  Indexed Access of Scalars:                ");
 363       a = fixed(is,ds,1);
 364       a(2) = fixed(is,ds,2);
 365       if (!isvector(a) || !isfixed(a))
 366         error("FAILED");
 367       endif
 368       a = fixed(is,ds,1);
 369       a(2) = fixed(is,ds,1i);
 370       if (!isvector(a) || !isfixed(a) || !iscomplex(a))
 371         error("FAILED");
 372       endif
 373       a = fixed(is,ds,1);
 374       a(1).dec = a(1).dec + 1;
 375       if (a.dec != ds+1)
 376         error("FAILED");
 377       endif
 378       a = fixed(is,ds,1i);
 379       a(2) = fixed(is,ds,2*1i);
 380       if (!isvector(a) || !isfixed(a) || !iscomplex(a))
 381         error("FAILED");
 382       endif
 383       a = fixed(is,ds,1i);
 384       a(1).dec = a(1).dec + 1;
 385       if (a.dec != ds*(1+1i)+1)
 386         error("FAILED");
 387       endif
 388       fprintf("PASSED\n");
 389       fprintf("  Indexed Access of Matrices:               ");
 390       x = (vals(2)-vals(1))*rand(n,n)+vals(1);
 391       a = fixed(is, ds, x);
 392       if (any(any(abs(a(1,:).x - x(1,:)) > feps)) || !isfixed(a))
 393         error ("FAILED");
 394       endif
 395       a(1,:).dec = a(1,:).dec + 1;
 396       if (a(1,:).dec != ds+1)
 397         error("FAILED");
 398       endif
 399       if (a(2,:).dec != ds)
 400         error("FAILED");
 401       endif
 402       a(1,1) = fixed(is,ds,1i);
 403       if (!ismatrix(a) || !isfixed(a) || !iscomplex(a))
 404         error("FAILED");
 405       endif
 406       x = ((vals(2)-vals(1))*rand(n,n)+vals(1) +
 407            ((vals(2)-vals(1))*rand(n,n)+vals(1)) * 1i);
 408       a = fixed(is, ds, x);
 409       if (any(any(abs(real(a(1,:).x) - real(x(1,:))) > feps)) ||
 410           any(any(abs(imag(a(1,:).x) - imag(x(1,:))) > feps)) ||
 411           !isfixed(a))
 412         error ("FAILED");
 413       endif
 414       a(1,:).dec = a(1,:).dec + 1;
 415       if (a(1,:).dec != ds*(1+1i)+1)
 416         error("FAILED");
 417       endif
 418       if (a(2,:).dec != ds*(1+1i))
 419         error("FAILED");
 420       endif
 421       fprintf("PASSED\n");
 422
 423       fprintf("\n<< Fixed Point Operators >>\n");
 424
 425       fprintf("  Logical Operators:                        ");
 426       if (_fixedpoint_test_bool_operator("==",is,ds,n) ||
 427           _fixedpoint_test_bool_operator("!=",is,ds,n) ||
 428           _fixedpoint_test_bool_operator(">",is,ds,n) ||
 429           _fixedpoint_test_bool_operator("<",is,ds,n) ||
 430           _fixedpoint_test_bool_operator(">=",is,ds,n) ||
 431           _fixedpoint_test_bool_operator("<=",is,ds,n))
 432         error("FAILED");
 433       else
 434         fprintf("PASSED\n");
 435       endif
 436
 437       fprintf("  Unary Operators:                          ");
 438       if (_fixedpoint_test_unary_operator("-",is,ds,n,1) ||
 439           _fixedpoint_test_unary_operator("+",is,ds,n,1) ||
 440           _fixedpoint_test_unary_operator("'",is,ds,n,0) ||
 441           _fixedpoint_test_unary_operator(".'",is,ds,n,0))
 442         error("FAILED");
 443       else
 444         fprintf("PASSED\n");
 445       endif
 446       # This has special meaning for fixed point. Whats the test !!!!
 447       # Is the damn operator even correct !!
 448       # b = ! a;
 449
 450       fprintf("  Arithmetic Operators:                     ");
 451       if (_fixedpoint_test_operator("-",is,ds,n,1,1,1) ||
 452           _fixedpoint_test_operator("+",is,ds,n,1,1,1) ||
 453           _fixedpoint_test_operator(".*",is,ds,n,1,1,1) ||
 454           _fixedpoint_test_operator("*",is,ds,n,1,1,1) ||
 455           _fixedpoint_test_operator("./",is,ds,n,0,1,1) ||
 456
 457 ##  Bug in octave 2.1.50 and earlier for el_ldiv in op-cs-s.cc.
 458           ((vers(1) > 2 || (vers(1) >= 2 && vers(2) > 1) ||
 459           (vers(1) >= 2 && vers(2) >= 1 && vers(3) > 50)) &&
 460            _fixedpoint_test_operator(".\\",is,ds,n,0,1,1)) ||
 461
 462 ## Rounding errors in complex pow functions make these fail
 463 ##        _fixedpoint_test_operator(".**",is,ds,n,0,1,1) ||
 464 ##        _fixedpoint_test_operator(".^",is,ds,n,0,1,1)
 465
 466 ## Can't do "matrix-by-matrix", as well as above problem
 467 ##        _fixedpoint_test_operator("**",is,ds,n,0,1,0) ||
 468 ##        _fixedpoint_test_operator("^",is,ds,n,0,1,0)  ||
 469
 470 ## Can't divide by a matrix as don't have inversion
 471           _fixedpoint_test_operator("/",is,ds,n,0,1,0) ||
 472           _fixedpoint_test_operator("\\",is,ds,n,0,0,1))
 473         error("FAILED");
 474       else
 475         fprintf("PASSED\n");
 476       endif
 477
 478       fprintf("\n<< Fixed Point Functions >>\n");
 479       fprintf("  Rounding Functions:                       ");
 480
 481       if (_fixedpoint_test_function("ceil",is,ds,n,0) ||
 482           _fixedpoint_test_function("floor",is,ds,n,0) ||
 483           _fixedpoint_test_function("round",is,ds,n,0))
 484         error("FAILED");
 485       else
 486         fprintf("PASSED\n");
 487       endif
 488
 489       fprintf("  Matrix Sum and Product Functions:         ");
 490       if (_fixedpoint_test_function("sum",is,ds,n,0) ||
 491           _fixedpoint_test_function("cumsum",is,ds,n,0))
 492         error("FAILED");
 493       endif
 494
 495       ## These can easily under- or  over-flow. Trick the code by
 496       ## foring ds=0, which simplifies the internal calculations,
 497       ## effectively limiting the values used to 1 and 0. If you
 498       ## can think of a better way, tell me...
 499       tmpds = ds;
 500       ds = 0;
 501       if (_fixedpoint_test_function("sumsq",is,ds,n,feps) ||
 502           _fixedpoint_test_function("prod",is,ds,n,feps) ||
 503           _fixedpoint_test_function("cumprod",is,ds,n,feps))
 504         error("FAILED");
 505       else
 506         fprintf("PASSED\n");
 507       endif
 508       ds = tmpds;
 509
 510       fprintf("  Miscellaneous Functions:                  ");
 511       if (_fixedpoint_test_function("real",is,ds,n,0) ||
 512           _fixedpoint_test_function("imag",is,ds,n,0) ||
 513           _fixedpoint_test_function("conj",is,ds,n,0) ||
 514           _fixedpoint_test_function("arg",is,ds,n,feps) ||
 515           _fixedpoint_test_function("angle",is,ds,n,feps) ||
 516           _fixedpoint_test_function("abs",is,ds,n,sqrt(2)*feps) ||
 517           _fixedpoint_test_function("sqrt",is,ds,n,sqrt(2)*feps))
 518         error("FAILED");
 519       else
 520         fprintf("PASSED\n");
 521       endif
 522
 523       ## These will underflow/overflow, use larger feps that should be
 524       ## ok upto ds = 6. Again if you can think of a better way to check!!
 525       fprintf("  Exponential Functions:                    ");
 526       if (_fixedpoint_test_function("log",is,ds,n,2*feps,0) ||
 527           _fixedpoint_test_function("log10",is,ds,n,20*feps,0) ||
 528           _fixedpoint_test_function("exp",is,ds,n,4*feps))
 529         error("FAILED");
 530       else
 531         fprintf("PASSED\n");
 532       endif
 533
 534       ## These will underflow/overflow, use larger feps that should be
 535       ## ok upto ds = 6. Again if you can think of a better way to check!!
 536       fprintf("  Trigonometric Functions:                  ");
 537       if (_fixedpoint_test_function("sin",is,ds,n,4*feps) ||
 538           _fixedpoint_test_function("sinh",is,ds,n,4*feps) ||
 539           _fixedpoint_test_function("cos",is,ds,n,4*feps) ||
 540           _fixedpoint_test_function("cosh",is,ds,n,4*feps) ||
 541           _fixedpoint_test_function("tan",is,ds,n,4*feps) ||
 542           _fixedpoint_test_function("tanh",is,ds,n,4*feps))
 543         error("FAILED");
 544       endif
 545
 546       ## Special case for fatan2
 547       fzero = fixed(is,ds,0);
 548       fzeros = fixed(is,ds,zeros(n,1));
 549       fone = fixed(is,ds,1);
 550       fones = fixed(is,ds,ones(n,1));
 551       fpi2 = fixed(is,ds,pi/2);
 552       fpi4 = fixed(is,ds,pi/4);
 553       if (fatan2(fzero,fzero) != fzero ||
 554           fatan2(fzero,fone) != fzero ||
 555           fatan2(fone,fzero) != fpi2 ||
 556           fatan2(fone,fone) != fpi4 ||
 557           any(fatan2(fzeros,fzeros) != fzero) ||
 558           any(fatan2(fzeros,fones) != fzero) ||
 559           any(fatan2(fones,fzeros) != fpi2) ||
 560           any(fatan2(fones,fones) != fpi4))
 561         error("FAILED");
 562       else
 563         fprintf("PASSED\n");
 564       endif
 565
 566       fprintf("\n");
 567     unwind_protect_cleanup
 568       page_screen_output(pso);
 569     end_unwind_protect
 570   else
 571     usage("fixedpoint: Unknown argument");
 572   endif
 573 endfunction
 574
 575 function ret = _fixedpoint_test_operator(op,is,ds,n,includezero,matrixleft,
 576                                          matrixright)
 577   # Properly test all 16 possible combinations of an operator
 578
 579   ret = 0;
 580
 581   tzero = 0;
 582   tone = 1;
 583   tzeros = zeros(n,n);
 584   tones = ones(n,n);
 585   tczero = 1i;
 586   tczeros = 1i*ones(n,n);
 587   tcone = 1 + 1i;
 588   tcones = (1 + 1i)*ones(n,n);
 589
 590   fzero = fixed(is,ds,tzero);
 591   fone = fixed(is,ds,tone);
 592   fzeros = fixed(is,ds,tzeros);
 593   fones = fixed(is,ds,tones);
 594   fczero = fixed(is,ds,tczero);
 595   fczeros = fixed(is,ds,tczeros);
 596   fcone = fixed(is,ds,tcone);
 597   fcones = fixed(is,ds,tcones);
 598
 599   # scalar by scalar
 600   if (includezero)
 601     t1 = eval(["tzero" op "tzero;"]);
 602     f1 = eval(["fzero" op "fzero;"]);
 603     if (t1 != f1.x)
 604       ret = 1;
 605       return;
 606     endif
 607     t1 = eval(["tzero" op "tone;"]);
 608     f1 = eval(["fzero" op "fone;"]);
 609     if (t1 != f1.x)
 610       ret = 1;
 611       return;
 612     endif
 613     t1 = eval(["tone" op "tzero;"]);
 614     f1 = eval(["fone" op "fzero;"]);
 615     if (t1 != f1.x)
 616       ret = 1;
 617       return;
 618     endif
 619     endif
 620   t1 = eval(["tone" op "tone;"]);
 621   f1 = eval(["fone" op "fone;"]);
 622   if (t1 != f1.x)
 623     ret = 1;
 624     return;
 625   endif
 626
 627   # scalar by matrix
 628   if (matrixright)
 629     if (includezero)
 630       t1 = eval(["tzero" op "tzeros;"]);
 631       f1 = eval(["fzero" op "fzeros;"]);
 632       if any(any(t1 != f1.x))
 633         ret = 1;
 634         return;
 635       endif
 636       t1 = eval(["tzero" op "tones;"]);
 637       f1 = eval(["fzero" op "fones;"]);
 638       if any(any(t1 != f1.x))
 639         ret = 1;
 640         return;
 641       endif
 642       t1 = eval(["tone" op "tzeros;"]);
 643       f1 = eval(["fone" op "fzeros;"]);
 644       if any(any(t1 != f1.x))
 645         ret = 1;
 646         return;
 647       endif
 648     endif
 649     t1 = eval(["tone" op "tones;"]);
 650     f1 = eval(["fone" op "fones;"]);
 651     if any(any(t1 != f1.x))
 652       ret = 1;
 653       return;
 654     endif
 655   endif
 656
 657   # matrix by scalar
 658   if (matrixleft)
 659     if (includezero)
 660       t1 = eval(["tzeros" op "tzero;"]);
 661       f1 = eval(["fzeros" op "fzero;"]);
 662       if any(any(t1 != f1.x))
 663         ret = 1;
 664         return;
 665       endif
 666       t1 = eval(["tzeros" op "tone;"]);
 667       f1 = eval(["fzeros" op "fone;"]);
 668       if any(any(t1 != f1.x))
 669         ret = 1;
 670         return;
 671       endif
 672       t1 = eval(["tones" op "tzero;"]);
 673       f1 = eval(["fones" op "fzero;"]);
 674       if any(any(t1 != f1.x))
 675         ret = 1;
 676         return;
 677       endif
 678     endif
 679     t1 = eval(["tones" op "tone;"]);
 680     f1 = eval(["fones" op "fone;"]);
 681     if any(any(t1 != f1.x))
 682       ret = 1;
 683       return;
 684     endif
 685   endif
 686
 687   # matrix by matrix
 688   if (matrixleft && matrixright)
 689     if (includezero)
 690       t1 = eval(["tzeros" op "tzeros;"]);
 691       f1 = eval(["fzeros" op "fzeros;"]);
 692       if any(any(t1 != f1.x))
 693         ret = 1;
 694         return;
 695       endif
 696       t1 = eval(["tzeros" op "tones;"]);
 697       f1 = eval(["fzeros" op "fones;"]);
 698       if any(any(t1 != f1.x))
 699         ret = 1;
 700         return;
 701       endif
 702       t1 = eval(["tones" op "tzeros;"]);
 703       f1 = eval(["fones" op "fzeros;"]);
 704       if any(any(t1 != f1.x))
 705         ret = 1;
 706         return;
 707       endif
 708     endif
 709     t1 = eval(["tones" op "tones;"]);
 710     f1 = eval(["fones" op "fones;"]);
 711     if any(any(t1 != f1.x))
 712       ret = 1;
 713       return;
 714     endif
 715   endif
 716
 717   # scalar by complex scalar
 718   if (includezero)
 719     t1 = eval(["tzero" op "tczero;"]);
 720     f1 = eval(["fzero" op "fczero;"]);
 721     if (t1 != f1.x)
 722       ret = 1;
 723       return;
 724     endif
 725     t1 = eval(["tzero" op "tcone;"]);
 726     f1 = eval(["fzero" op "fcone;"]);
 727     if (t1 != f1.x)
 728       ret = 1;
 729       return;
 730     endif
 731   endif
 732   t1 = eval(["tone" op "tczero;"]);
 733   f1 = eval(["fone" op "fczero;"]);
 734   if (t1 != f1.x)
 735     ret = 1;
 736     return;
 737   endif
 738   t1 = eval(["tone" op "tcone;"]);
 739   f1 = eval(["fone" op "fcone;"]);
 740   if (t1 != f1.x)
 741     ret = 1;
 742     return;
 743   endif
 744
 745   # scalar by complex matrix
 746   if (matrixright)
 747     if (includezero)
 748       t1 = eval(["tzero" op "tczeros;"]);
 749       f1 = eval(["fzero" op "fczeros;"]);
 750       if any(any(t1 != f1.x))
 751         ret = 1;
 752         return;
 753       endif
 754       t1 = eval(["tzero" op "tcones;"]);
 755       f1 = eval(["fzero" op "fcones;"]);
 756       if any(any(t1 != f1.x))
 757         ret = 1;
 758         return;
 759       endif
 760     endif
 761     t1 = eval(["tone" op "tczeros;"]);
 762     f1 = eval(["fone" op "fczeros;"]);
 763     if any(any(t1 != f1.x))
 764       ret = 1;
 765       return;
 766     endif
 767     t1 = eval(["tone" op "tcones;"]);
 768     f1 = eval(["fone" op "fcones;"]);
 769     if any(any(t1 != f1.x))
 770       ret = 1;
 771       return;
 772     endif
 773   endif
 774
 775   # matrix by complex scalar
 776   if (matrixleft)
 777     if (includezero)
 778       t1 = eval(["tzeros" op "tczero;"]);
 779       f1 = eval(["fzeros" op "fczero;"]);
 780       if any(any(t1 != f1.x))
 781         ret = 1;
 782         return;
 783       endif
 784       t1 = eval(["tzeros" op "tcone;"]);
 785       f1 = eval(["fzeros" op "fcone;"]);
 786       if any(any(t1 != f1.x))
 787         ret = 1;
 788         return;
 789       endif
 790     endif
 791     t1 = eval(["tones" op "tczero;"]);
 792     f1 = eval(["fones" op "fczero;"]);
 793     if any(any(t1 != f1.x))
 794       ret = 1;
 795       return;
 796     endif
 797     t1 = eval(["tones" op "tcone;"]);
 798     f1 = eval(["fones" op "fcone;"]);
 799     if any(any(t1 != f1.x))
 800       ret = 1;
 801       return;
 802     endif
 803   endif
 804
 805   # matrix by complex matrix
 806   if (matrixleft && matrixright)
 807     if (includezero)
 808       t1 = eval(["tzeros" op "tczeros;"]);
 809       f1 = eval(["fzeros" op "fczeros;"]);
 810       if any(any(t1 != f1.x))
 811         ret = 1;
 812         return;
 813       endif
 814       t1 = eval(["tzeros" op "tcones;"]);
 815       f1 = eval(["fzeros" op "fcones;"]);
 816       if any(any(t1 != f1.x))
 817         ret = 1;
 818         return;
 819       endif
 820     endif
 821     t1 = eval(["tones" op "tczeros;"]);
 822     f1 = eval(["fones" op "fczeros;"]);
 823     if any(any(t1 != f1.x))
 824       ret = 1;
 825       return;
 826     endif
 827     t1 = eval(["tones" op "tcones;"]);
 828     f1 = eval(["fones" op "fcones;"]);
 829     if any(any(t1 != f1.x))
 830       ret = 1;
 831       return;
 832     endif
 833   endif
 834
 835   # complex scalar by scalar
 836   if (includezero)
 837     t1 = eval(["tczero" op "tzero;"]);
 838     f1 = eval(["fczero" op "fzero;"]);
 839     if (t1 != f1.x)
 840       ret = 1;
 841       return;
 842     endif
 843     t1 = eval(["tcone" op "tzero;"]);
 844     f1 = eval(["fcone" op "fzero;"]);
 845     if (t1 != f1.x)
 846       ret = 1;
 847       return;
 848     endif
 849   endif
 850   t1 = eval(["tczero" op "tone;"]);
 851   f1 = eval(["fczero" op "fone;"]);
 852   if (t1 != f1.x)
 853     ret = 1;
 854     return;
 855   endif
 856   t1 = eval(["tcone" op "tone;"]);
 857   f1 = eval(["fcone" op "fone;"]);
 858   if (t1 != f1.x)
 859     ret = 1;
 860     return;
 861   endif
 862
 863   # complex scalar by matrix
 864   if (matrixright)
 865     if (includezero)
 866       t1 = eval(["tczero" op "tzeros;"]);
 867       f1 = eval(["fczero" op "fzeros;"]);
 868       if any(any(t1 != f1.x))
 869         ret = 1;
 870         return;
 871       endif
 872       t1 = eval(["tcone" op "tzeros;"]);
 873       f1 = eval(["fcone" op "fzeros;"]);
 874       if any(any(t1 != f1.x))
 875         ret = 1;
 876         return;
 877       endif
 878     endif
 879     t1 = eval(["tczero" op "tones;"]);
 880     f1 = eval(["fczero" op "fones;"]);
 881     if any(any(t1 != f1.x))
 882       ret = 1;
 883       return;
 884     endif
 885     t1 = eval(["tcone" op "tones;"]);
 886     f1 = eval(["fcone" op "fones;"]);
 887     if any(any(t1 != f1.x))
 888       ret = 1;
 889       return;
 890     endif
 891   endif
 892
 893   # complex matrix by scalar
 894   if (matrixleft)
 895     if (includezero)
 896       t1 = eval(["tczeros" op "tzero;"]);
 897       f1 = eval(["fczeros" op "fzero;"]);
 898       if any(any(t1 != f1.x))
 899         ret = 1;
 900         return;
 901       endif
 902       t1 = eval(["tcones" op "tzero;"]);
 903       f1 = eval(["fcones" op "fzero;"]);
 904       if any(any(t1 != f1.x))
 905         ret = 1;
 906         return;
 907       endif
 908     endif
 909     t1 = eval(["tczeros" op "tone;"]);
 910     f1 = eval(["fczeros" op "fone;"]);
 911     if any(any(t1 != f1.x))
 912       ret = 1;
 913       return;
 914     endif
 915     t1 = eval(["tcones" op "tone;"]);
 916     f1 = eval(["fcones" op "fone;"]);
 917     if any(any(t1 != f1.x))
 918       ret = 1;
 919       return;
 920     endif
 921   endif
 922
 923   # complex matrix by matrix
 924   if (matrixleft && matrixright)
 925     if (includezero)
 926       t1 = eval(["tczeros" op "tzeros;"]);
 927       f1 = eval(["fczeros" op "fzeros;"]);
 928       if any(any(t1 != f1.x))
 929         ret = 1;
 930         return;
 931       endif
 932       t1 = eval(["tcones" op "tzeros;"]);
 933       f1 = eval(["fcones" op "fzeros;"]);
 934       if any(any(t1 != f1.x))
 935         ret = 1;
 936         return;
 937       endif
 938     endif
 939     t1 = eval(["tczeros" op "tones;"]);
 940     f1 = eval(["fczeros" op "fones;"]);
 941     if any(any(t1 != f1.x))
 942       ret = 1;
 943       return;
 944     endif
 945     t1 = eval(["tcones" op "tones;"]);
 946     f1 = eval(["fcones" op "fones;"]);
 947     if any(any(t1 != f1.x))
 948       ret = 1;
 949       return;
 950     endif
 951   endif
 952
 953   # complex scalar by complex scalar
 954   t1 = eval(["tczero" op "tczero;"]);
 955   f1 = eval(["fczero" op "fczero;"]);
 956   if (t1 != f1.x)
 957     ret = 1;
 958     return;
 959   endif
 960   t1 = eval(["tczero" op "tcone;"]);
 961   f1 = eval(["fczero" op "fcone;"]);
 962   if (t1 != f1.x)
 963     ret = 1;
 964     return;
 965   endif
 966   t1 = eval(["tcone" op "tczero;"]);
 967   f1 = eval(["fcone" op "fczero;"]);
 968   if (t1 != f1.x)
 969     ret = 1;
 970     return;
 971   endif
 972   t1 = eval(["tcone" op "tcone;"]);
 973   f1 = eval(["fcone" op "fcone;"]);
 974   if (t1 != f1.x)
 975     ret = 1;
 976     return;
 977   endif
 978
 979   # complex scalar by complex matrix
 980   if (matrixright)
 981     t1 = eval(["tczero" op "tczeros;"]);
 982     f1 = eval(["fczero" op "fczeros;"]);
 983     if any(any(t1 != f1.x))
 984       ret = 1;
 985       return;
 986     endif
 987     t1 = eval(["tczero" op "tcones;"]);
 988     f1 = eval(["fczero" op "fcones;"]);
 989     if any(any(t1 != f1.x))
 990       ret = 1;
 991       return;
 992     endif
 993     t1 = eval(["tcone" op "tczeros;"]);
 994     f1 = eval(["fcone" op "fczeros;"]);
 995     if any(any(t1 != f1.x))
 996       ret = 1;
 997       return;
 998     endif
 999     t1 = eval(["tcone" op "tcones;"]);
1000     f1 = eval(["fcone" op "fcones;"]);
1001     if any(any(t1 != f1.x))
1002       ret = 1;
1003       return;
1004     endif
1005   endif
1006
1007   # complex matrix by complex scalar
1008   if (matrixleft)
1009     t1 = eval(["tczeros" op "tczero;"]);
1010     f1 = eval(["fczeros" op "fczero;"]);
1011     if any(any(t1 != f1.x))
1012       ret = 1;
1013       return;
1014     endif
1015     t1 = eval(["tczeros" op "tcone;"]);
1016     f1 = eval(["fczeros" op "fcone;"]);
1017     if any(any(t1 != f1.x))
1018       ret = 1;
1019       return;
1020     endif
1021     t1 = eval(["tcones" op "tczero;"]);
1022     f1 = eval(["fcones" op "fczero;"]);
1023     if any(any(t1 != f1.x))
1024       ret = 1;
1025       return;
1026     endif
1027     t1 = eval(["tcones" op "tcone;"]);
1028     f1 = eval(["fcones" op "fcone;"]);
1029     if any(any(t1 != f1.x))
1030       ret = 1;
1031       return;
1032     endif
1033   endif
1034
1035   # complex matrix by complex matrix
1036   if (matrixleft && matrixright)
1037     t1 = eval(["tczeros" op "tczeros;"]);
1038     f1 = eval(["fczeros" op "fczeros;"]);
1039     if any(any(t1 != f1.x))
1040       ret = 1;
1041       return;
1042     endif
1043     t1 = eval(["tczeros" op "tcones;"]);
1044     f1 = eval(["fczeros" op "fcones;"]);
1045     if any(any(t1 != f1.x))
1046       ret = 1;
1047       return;
1048     endif
1049     t1 = eval(["tcones" op "tczeros;"]);
1050     f1 = eval(["fcones" op "fczeros;"]);
1051     if any(any(t1 != f1.x))
1052       ret = 1;
1053       return;
1054     endif
1055     t1 = eval(["tcones" op "tcones;"]);
1056     f1 = eval(["fcones" op "fcones;"]);
1057     if any(any(t1 != f1.x))
1058       ret = 1;
1059       return;
1060     endif
1061   endif
1062
1063 endfunction
1064
1065 function ret = _fixedpoint_test_bool_operator(op,is,ds,n)
1066   # Properly test all 16 possible combinations of an operator
1067
1068   ret = 0;
1069
1070   tzero = 0;
1071   tone = 1;
1072   tzeros = zeros(n,n);
1073   tones = ones(n,n);
1074   tczero = 1i;
1075   tczeros = 1i*ones(n,n);
1076   tcone = 1 + 1i;
1077   tcones = (1 + 1i)*ones(n,n);
1078
1079   fzero = fixed(is,ds,tzero);
1080   fone = fixed(is,ds,tone);
1081   fzeros = fixed(is,ds,tzeros);
1082   fones = fixed(is,ds,tones);
1083   fczero = fixed(is,ds,tczero);
1084   fczeros = fixed(is,ds,tczeros);
1085   fcone = fixed(is,ds,tcone);
1086   fcones = fixed(is,ds,tcones);
1087
1088   # scalar by scalar
1089   t1 = eval(["tzero" op "tzero;"]);
1090   f1 = eval(["fzero" op "fzero;"]);
1091   if (t1 != f1)
1092     ret = 1;
1093     return;
1094   endif
1095   t1 = eval(["tzero" op "tone;"]);
1096   f1 = eval(["fzero" op "fone;"]);
1097   if (t1 != f1)
1098     ret = 1;
1099     return;
1100   endif
1101   t1 = eval(["tone" op "tzero;"]);
1102   f1 = eval(["fone" op "fzero;"]);
1103   if (t1 != f1)
1104     ret = 1;
1105     return;
1106   endif
1107   t1 = eval(["tone" op "tone;"]);
1108   f1 = eval(["fone" op "fone;"]);
1109   if (t1 != f1)
1110     ret = 1;
1111     return;
1112   endif
1113
1114   # scalar by matrix
1115   t1 = eval(["tzero" op "tzeros;"]);
1116   f1 = eval(["fzero" op "fzeros;"]);
1117   if any(any(t1 != f1))
1118     ret = 1;
1119     return;
1120   endif
1121   t1 = eval(["tzero" op "tones;"]);
1122   f1 = eval(["fzero" op "fones;"]);
1123   if any(any(t1 != f1))
1124     ret = 1;
1125     return;
1126   endif
1127   t1 = eval(["tone" op "tzeros;"]);
1128   f1 = eval(["fone" op "fzeros;"]);
1129   if any(any(t1 != f1))
1130     ret = 1;
1131     return;
1132   endif
1133   t1 = eval(["tone" op "tones;"]);
1134   f1 = eval(["fone" op "fones;"]);
1135   if any(any(t1 != f1))
1136     ret = 1;
1137     return;
1138   endif
1139
1140   # matrix by scalar
1141   t1 = eval(["tzeros" op "tzero;"]);
1142   f1 = eval(["fzeros" op "fzero;"]);
1143   if any(any(t1 != f1))
1144     ret = 1;
1145     return;
1146   endif
1147   t1 = eval(["tzeros" op "tone;"]);
1148   f1 = eval(["fzeros" op "fone;"]);
1149   if any(any(t1 != f1))
1150     ret = 1;
1151     return;
1152   endif
1153   t1 = eval(["tones" op "tzero;"]);
1154   f1 = eval(["fones" op "fzero;"]);
1155   if any(any(t1 != f1))
1156     ret = 1;
1157     return;
1158   endif
1159   t1 = eval(["tones" op "tone;"]);
1160   f1 = eval(["fones" op "fone;"]);
1161   if any(any(t1 != f1))
1162     ret = 1;
1163     return;
1164   endif
1165
1166   # matrix by matrix
1167   t1 = eval(["tzeros" op "tzeros;"]);
1168   f1 = eval(["fzeros" op "fzeros;"]);
1169   if any(any(t1 != f1))
1170     ret = 1;
1171     return;
1172   endif
1173   t1 = eval(["tzeros" op "tones;"]);
1174   f1 = eval(["fzeros" op "fones;"]);
1175   if any(any(t1 != f1))
1176     ret = 1;
1177     return;
1178   endif
1179   t1 = eval(["tones" op "tzeros;"]);
1180   f1 = eval(["fones" op "fzeros;"]);
1181   if any(any(t1 != f1))
1182     ret = 1;
1183     return;
1184   endif
1185   t1 = eval(["tones" op "tones;"]);
1186   f1 = eval(["fones" op "fones;"]);
1187   if any(any(t1 != f1))
1188     ret = 1;
1189     return;
1190   endif
1191
1192
1193   # scalar by complex scalar
1194   t1 = eval(["tzero" op "tczero;"]);
1195   f1 = eval(["fzero" op "fczero;"]);
1196   if (t1 != f1)
1197     ret = 1;
1198     return;
1199   endif
1200   t1 = eval(["tzero" op "tcone;"]);
1201   f1 = eval(["fzero" op "fcone;"]);
1202   if (t1 != f1)
1203     ret = 1;
1204     return;
1205   endif
1206   t1 = eval(["tone" op "tczero;"]);
1207   f1 = eval(["fone" op "fczero;"]);
1208   if (t1 != f1)
1209     ret = 1;
1210     return;
1211   endif
1212   t1 = eval(["tone" op "tcone;"]);
1213   f1 = eval(["fone" op "fcone;"]);
1214   if (t1 != f1)
1215     ret = 1;
1216     return;
1217   endif
1218
1219   # scalar by complex matrix
1220   t1 = eval(["tzero" op "tczeros;"]);
1221   f1 = eval(["fzero" op "fczeros;"]);
1222   if any(any(t1 != f1))
1223     ret = 1;
1224     return;
1225   endif
1226   t1 = eval(["tzero" op "tcones;"]);
1227   f1 = eval(["fzero" op "fcones;"]);
1228   if any(any(t1 != f1))
1229     ret = 1;
1230     return;
1231   endif
1232   t1 = eval(["tone" op "tczeros;"]);
1233   f1 = eval(["fone" op "fczeros;"]);
1234   if any(any(t1 != f1))
1235     ret = 1;
1236     return;
1237   endif
1238   t1 = eval(["tone" op "tcones;"]);
1239   f1 = eval(["fone" op "fcones;"]);
1240   if any(any(t1 != f1))
1241     ret = 1;
1242     return;
1243   endif
1244
1245   # matrix by complex scalar
1246   t1 = eval(["tzeros" op "tczero;"]);
1247   f1 = eval(["fzeros" op "fczero;"]);
1248   if any(any(t1 != f1))
1249     ret = 1;
1250     return;
1251   endif
1252   t1 = eval(["tzeros" op "tcone;"]);
1253   f1 = eval(["fzeros" op "fcone;"]);
1254   if any(any(t1 != f1))
1255     ret = 1;
1256     return;
1257   endif
1258   t1 = eval(["tones" op "tczero;"]);
1259   f1 = eval(["fones" op "fczero;"]);
1260   if any(any(t1 != f1))
1261     ret = 1;
1262     return;
1263   endif
1264   t1 = eval(["tones" op "tcone;"]);
1265   f1 = eval(["fones" op "fcone;"]);
1266   if any(any(t1 != f1))
1267     ret = 1;
1268     return;
1269   endif
1270
1271   # matrix by complex matrix
1272   t1 = eval(["tzeros" op "tczeros;"]);
1273   f1 = eval(["fzeros" op "fczeros;"]);
1274   if any(any(t1 != f1))
1275     ret = 1;
1276     return;
1277   endif
1278   t1 = eval(["tzeros" op "tcones;"]);
1279   f1 = eval(["fzeros" op "fcones;"]);
1280   if any(any(t1 != f1))
1281     ret = 1;
1282     return;
1283   endif
1284   t1 = eval(["tones" op "tczeros;"]);
1285   f1 = eval(["fones" op "fczeros;"]);
1286   if any(any(t1 != f1))
1287     ret = 1;
1288     return;
1289   endif
1290   t1 = eval(["tones" op "tcones;"]);
1291   f1 = eval(["fones" op "fcones;"]);
1292   if any(any(t1 != f1))
1293     ret = 1;
1294     return;
1295   endif
1296
1297   # complex scalar by scalar
1298   t1 = eval(["tczero" op "tzero;"]);
1299   f1 = eval(["fczero" op "fzero;"]);
1300   if (t1 != f1)
1301     ret = 1;
1302     return;
1303   endif
1304   t1 = eval(["tczero" op "tone;"]);
1305   f1 = eval(["fczero" op "fone;"]);
1306   if (t1 != f1)
1307     ret = 1;
1308     return;
1309   endif
1310   t1 = eval(["tcone" op "tzero;"]);
1311   f1 = eval(["fcone" op "fzero;"]);
1312   if (t1 != f1)
1313     ret = 1;
1314     return;
1315   endif
1316   t1 = eval(["tcone" op "tone;"]);
1317   f1 = eval(["fcone" op "fone;"]);
1318   if (t1 != f1)
1319     ret = 1;
1320     return;
1321   endif
1322
1323   # complex scalar by matrix
1324   t1 = eval(["tczero" op "tzeros;"]);
1325   f1 = eval(["fczero" op "fzeros;"]);
1326   if any(any(t1 != f1))
1327     ret = 1;
1328     return;
1329   endif
1330   t1 = eval(["tczero" op "tones;"]);
1331   f1 = eval(["fczero" op "fones;"]);
1332   if any(any(t1 != f1))
1333     ret = 1;
1334     return;
1335   endif
1336   t1 = eval(["tcone" op "tzeros;"]);
1337   f1 = eval(["fcone" op "fzeros;"]);
1338   if any(any(t1 != f1))
1339     ret = 1;
1340     return;
1341   endif
1342   t1 = eval(["tcone" op "tones;"]);
1343   f1 = eval(["fcone" op "fones;"]);
1344   if any(any(t1 != f1))
1345     ret = 1;
1346     return;
1347   endif
1348
1349   # complex matrix by scalar
1350   t1 = eval(["tczeros" op "tzero;"]);
1351   f1 = eval(["fczeros" op "fzero;"]);
1352   if any(any(t1 != f1))
1353     ret = 1;
1354     return;
1355   endif
1356   t1 = eval(["tczeros" op "tone;"]);
1357   f1 = eval(["fczeros" op "fone;"]);
1358   if any(any(t1 != f1))
1359     ret = 1;
1360     return;
1361   endif
1362   t1 = eval(["tcones" op "tzero;"]);
1363   f1 = eval(["fcones" op "fzero;"]);
1364   if any(any(t1 != f1))
1365     ret = 1;
1366     return;
1367   endif
1368   t1 = eval(["tcones" op "tone;"]);
1369   f1 = eval(["fcones" op "fone;"]);
1370   if any(any(t1 != f1))
1371     ret = 1;
1372     return;
1373   endif
1374
1375   # complex matrix by matrix
1376   t1 = eval(["tczeros" op "tzeros;"]);
1377   f1 = eval(["fczeros" op "fzeros;"]);
1378   if any(any(t1 != f1))
1379     ret = 1;
1380     return;
1381   endif
1382   t1 = eval(["tczeros" op "tones;"]);
1383   f1 = eval(["fczeros" op "fones;"]);
1384   if any(any(t1 != f1))
1385     ret = 1;
1386     return;
1387   endif
1388   t1 = eval(["tcones" op "tzeros;"]);
1389   f1 = eval(["fcones" op "fzeros;"]);
1390   if any(any(t1 != f1))
1391     ret = 1;
1392     return;
1393   endif
1394   t1 = eval(["tcones" op "tones;"]);
1395   f1 = eval(["fcones" op "fones;"]);
1396   if any(any(t1 != f1))
1397     ret = 1;
1398     return;
1399   endif
1400
1401   # complex scalar by complex scalar
1402   t1 = eval(["tczero" op "tczero;"]);
1403   f1 = eval(["fczero" op "fczero;"]);
1404   if (t1 != f1)
1405     ret = 1;
1406     return;
1407   endif
1408   t1 = eval(["tczero" op "tcone;"]);
1409   f1 = eval(["fczero" op "fcone;"]);
1410   if (t1 != f1)
1411     ret = 1;
1412     return;
1413   endif
1414   t1 = eval(["tcone" op "tczero;"]);
1415   f1 = eval(["fcone" op "fczero;"]);
1416   if (t1 != f1)
1417     ret = 1;
1418     return;
1419   endif
1420   t1 = eval(["tcone" op "tcone;"]);
1421   f1 = eval(["fcone" op "fcone;"]);
1422   if (t1 != f1)
1423     ret = 1;
1424     return;
1425   endif
1426
1427   # complex scalar by complex matrix
1428   t1 = eval(["tczero" op "tczeros;"]);
1429   f1 = eval(["fczero" op "fczeros;"]);
1430   if any(any(t1 != f1))
1431     ret = 1;
1432     return;
1433   endif
1434   t1 = eval(["tczero" op "tcones;"]);
1435   f1 = eval(["fczero" op "fcones;"]);
1436   if any(any(t1 != f1))
1437     ret = 1;
1438     return;
1439   endif
1440   t1 = eval(["tcone" op "tczeros;"]);
1441   f1 = eval(["fcone" op "fczeros;"]);
1442   if any(any(t1 != f1))
1443     ret = 1;
1444     return;
1445   endif
1446   t1 = eval(["tcone" op "tcones;"]);
1447   f1 = eval(["fcone" op "fcones;"]);
1448   if any(any(t1 != f1))
1449     ret = 1;
1450     return;
1451   endif
1452
1453   # complex matrix by complex scalar
1454   t1 = eval(["tczeros" op "tczero;"]);
1455   f1 = eval(["fczeros" op "fczero;"]);
1456   if any(any(t1 != f1))
1457     ret = 1;
1458     return;
1459   endif
1460   t1 = eval(["tczeros" op "tcone;"]);
1461   f1 = eval(["fczeros" op "fcone;"]);
1462   if any(any(t1 != f1))
1463     ret = 1;
1464     return;
1465   endif
1466   t1 = eval(["tcones" op "tczero;"]);
1467   f1 = eval(["fcones" op "fczero;"]);
1468   if any(any(t1 != f1))
1469     ret = 1;
1470     return;
1471   endif
1472   t1 = eval(["tcones" op "tcone;"]);
1473   f1 = eval(["fcones" op "fcone;"]);
1474   if any(any(t1 != f1))
1475     ret = 1;
1476     return;
1477   endif
1478
1479   # complex matrix by complex matrix
1480   t1 = eval(["tczeros" op "tczeros;"]);
1481   f1 = eval(["fczeros" op "fczeros;"]);
1482   if any(any(t1 != f1))
1483     ret = 1;
1484     return;
1485   endif
1486   t1 = eval(["tczeros" op "tcones;"]);
1487   f1 = eval(["fczeros" op "fcones;"]);
1488   if any(any(t1 != f1))
1489     ret = 1;
1490     return;
1491   endif
1492   t1 = eval(["tcones" op "tczeros;"]);
1493   f1 = eval(["fcones" op "fczeros;"]);
1494   if any(any(t1 != f1))
1495     ret = 1;
1496     return;
1497   endif
1498   t1 = eval(["tcones" op "tcones;"]);
1499   f1 = eval(["fcones" op "fcones;"]);
1500   if any(any(t1 != f1))
1501     ret = 1;
1502     return;
1503   endif
1504
1505 endfunction
1506
1507 function ret = _fixedpoint_test_unary_operator(op,is,ds,n,front)
1508   # Properly test all 4 possible combinations of an operator
1509
1510   ret = 0;
1511
1512   tzero = 0;
1513   tone = 1;
1514   tzeros = zeros(n,n);
1515   tones = ones(n,n);
1516   tczero = 1i;
1517   tczeros = 1i*ones(n,n);
1518   tcone = 1 + 1i;
1519   tcones = (1 + 1i)*ones(n,n);
1520
1521   fzero = fixed(is,ds,tzero);
1522   fone = fixed(is,ds,tone);
1523   fzeros = fixed(is,ds,tzeros);
1524   fones = fixed(is,ds,tones);
1525   fczero = fixed(is,ds,tczero);
1526   fczeros = fixed(is,ds,tczeros);
1527   fcone = fixed(is,ds,tcone);
1528   fcones = fixed(is,ds,tcones);
1529
1530   if (front)
1531     op1 = op;
1532     op2 = "";
1533   else
1534     op1 = "";
1535     op2 = op;
1536   endif
1537
1538   # scalar by scalar
1539   t1 = eval([ op1 "tzero" op2 ";"]);
1540   f1 = eval([ op1 "fzero" op2 ";"]);
1541   if (t1 != f1.x)
1542     ret = 1;
1543     return;
1544   endif
1545   t1 = eval([ op1 "tone" op2 ";"]);
1546   f1 = eval([ op1 "fone" op2 ";"]);
1547   if (t1 != f1.x)
1548     ret = 1;
1549     return;
1550   endif
1551
1552   # matrix
1553   t1 = eval([ op1 "tzeros" op2 ";"]);
1554   f1 = eval([ op1 "fzeros" op2 ";"]);
1555   if any(any(t1 != f1.x))
1556     ret = 1;
1557     return;
1558   endif
1559   t1 = eval([ op1 "tones" op2 ";"]);
1560   f1 = eval([ op1 "fones" op2 ";"]);
1561   if any(any(t1 != f1.x))
1562     ret = 1;
1563     return;
1564   endif
1565
1566   # complex scalar
1567   t1 = eval([ op1 "tczero" op2 ";"]);
1568   f1 = eval([ op1 "fczero" op2 ";"]);
1569   if (t1 != f1.x)
1570     ret = 1;
1571     return;
1572   endif
1573   t1 = eval([ op1 "tcone" op2 ";"]);
1574   f1 = eval([ op1 "fcone" op2 ";"]);
1575   if (t1 != f1.x)
1576     ret = 1;
1577     return;
1578   endif
1579
1580   # complex matrix
1581   t1 = eval([ op1 "tczeros" op2 ";"]);
1582   f1 = eval([ op1 "fczeros" op2 ";"]);
1583   if any(any(t1 != f1.x))
1584     ret = 1;
1585     return;
1586   endif
1587   t1 = eval([ op1 "tcones" op2 ";"]);
1588   f1 = eval([ op1 "fcones" op2 ";"]);
1589   if any(any(t1 != f1.x))
1590     ret = 1;
1591     return;
1592   endif
1593
1594 endfunction
1595
1596 function ret = _fixedpoint_test_function(func,is,ds,n,eps,includezero)
1597   # Properly test all 4 possible combinations of a function
1598
1599   if (nargin < 6)
1600     includezero = 1;
1601   endif
1602
1603   ret = 0;
1604   tzero = 0;
1605   tone = 1;
1606   tzeros = zeros(n,n);
1607   tones = ones(n,n);
1608   tczero = 1i;
1609   tczeros = 1i*ones(n,n);
1610   tcone = 1 + 1i;
1611   tcones = (1 + 1i)*ones(n,n);
1612   if (ds == 0)
1613     tval = tone;
1614     tvals = tones;
1615     tcval = tcone;
1616     tcvals = tcones;
1617   else
1618     tval = 0.5;
1619     if (includezero)
1620       tvals = ones(n,1) * (round(2^ds * (0:(n-1))/(n-1)) / 2^ds);
1621     else
1622       tvals = ones(n,1) * (round(2^ds * (1 + 2^(is-1)*(1:(n-1))/(n-1)))
1623                            / 2^ds);
1624     endif
1625     tcval = tval * (1+1i);
1626     tcvals = tvals * (1+1i);
1627   endif
1628
1629   fzero = fixed(is,ds,tzero);
1630   fone = fixed(is,ds,tone);
1631   fzeros = fixed(is,ds,tzeros);
1632   fones = fixed(is,ds,tones);
1633   fczero = fixed(is,ds,tczero);
1634   fczeros = fixed(is,ds,tczeros);
1635   fcone = fixed(is,ds,tcone);
1636   fcones = fixed(is,ds,tcones);
1637   fval = fixed(is,ds,tval);
1638   fvals = fixed(is,ds,tvals);
1639   fcval = fixed(is,ds,tcval);
1640   fcvals = fixed(is,ds,tcvals);
1641
1642   ffunc = [ "f" func];
1643
1644   # scalar by scalar
1645   if (includezero)
1646     t1 = eval([ func "(tzero);"]);
1647     f1 = eval([ ffunc "(fzero);"]);
1648     if (abs(t1 - f1.x) > eps)
1649       ret = 1;
1650       return;
1651     endif
1652   endif
1653   t1 = eval([ func "(tone);"]);
1654   f1 = eval([ ffunc "(fone);"]);
1655   if (abs(t1 - f1.x) > eps)
1656     ret = 1;
1657     return;
1658   endif
1659   t1 = eval([ func "(tval);"]);
1660   f1 = eval([ ffunc "(fval);"]);
1661   if (abs(t1 - f1.x) > eps)
1662     ret = 1;
1663     return;
1664   endif
1665
1666   # matrix
1667   if (includezero)
1668     t1 = eval([ func "(tzeros);"]);
1669     f1 = eval([ ffunc "(fzeros);"]);
1670     if any(any(abs(t1 - f1.x) > eps))
1671       ret = 1;
1672       return;
1673     endif
1674   endif
1675   t1 = eval([ func "(tones);"]);
1676   f1 = eval([ ffunc "(fones);"]);
1677   if any(any(abs(t1 - f1.x) > eps))
1678     ret = 1;
1679     return;
1680   endif
1681   t1 = eval([ func "(tvals);"]);
1682   f1 = eval([ ffunc "(fvals);"]);
1683   if any(any(abs(t1 - f1.x) > eps))
1684     ret = 1;
1685     return;
1686   endif
1687
1688   # complex scalar
1689   if (includezero)
1690     t1 = eval([ func "(tczero);"]);
1691     f1 = eval([ ffunc "(fczero);"]);
1692     if ((abs(real(t1) - freal(f1).x) > eps) ||
1693         (abs(imag(t1) - fimag(f1).x) > eps))
1694       ret = 1;
1695       return;
1696     endif
1697   endif
1698   t1 = eval([ func "(tcone);"]);
1699   f1 = eval([ ffunc "(fcone);"]);
1700   if ((abs(real(t1) - freal(f1).x) > eps) || (abs(imag(t1) - fimag(f1).x) > eps))
1701     ret = 1;
1702     return;
1703   endif
1704   t1 = eval([ func "(tcval);"]);
1705   f1 = eval([ ffunc "(fcval);"]);
1706   if ((abs(real(t1) - freal(f1).x) > eps) || (abs(imag(t1) - fimag(f1).x) > eps))
1707     ret = 1;
1708     return;
1709   endif
1710
1711   # complex matrix
1712   if (includezero)
1713     t1 = eval([ func "(tczeros);"]);
1714     f1 = eval([ ffunc "(fczeros);"]);
1715     if (any(any(abs(real(t1) - freal(f1).x) > eps)) ||
1716         any(any(abs(imag(t1) - fimag(f1).x) > eps)))
1717       ret = 1;
1718       return;
1719     endif
1720   endif
1721   t1 = eval([ func "(tcones);"]);
1722   f1 = eval([ ffunc "(fcones);"]);
1723   if (any(any(abs(real(t1) - freal(f1).x) > eps)) ||
1724       any(any(abs(imag(t1) - fimag(f1).x) > eps)))
1725     ret = 1;
1726     return;
1727   endif
1728   t1 = eval([ func "(tcvals);"]);
1729   f1 = eval([ ffunc "(fcvals);"]);
1730   if (any(any(abs(real(t1) - freal(f1).x) > eps)) ||
1731       any(any(abs(imag(t1) - fimag(f1).x) > eps)))
1732     ret = 1;
1733     return;
1734   endif
1735
1736 endfunction
1737
1738 %!test
1739 %! try fixedpoint("test");
1740 %! catch disp(lasterr()); end
1741
1742 ## Local Variables: ***
1743 ## mode: Octave ***
1744 ## End: ***