// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
-// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
-// <-- JVM NOT MANDATORY -->
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+// <-- CLI SHELL MODE -->
// <-- ENGLISH IMPOSED -->
function flag = MY_assert_equal ( computed , expected )
- if computed==expected then
- flag = 1;
- else
- flag = 0;
- end
- if flag <> 1 then bugmes();quit;end
+ if computed==expected then
+ flag = 1;
+ else
+ flag = 0;
+ end
+ if flag <> 1 then bugmes();quit;end
endfunction
format("v",10);
//
-dmax = -log10(2^(-53));
+dmax = -log(2^(-53))/log(10);
//
computed = assert_computedigits ( 1 , 1 );
MY_assert_equal ( computed , dmax );
MY_assert_equal ( computed , 0 );
//
computed = assert_computedigits ( 3.1415926 , %pi );
-MY_assert_equal ( computed , 7.76806779280040160 );
+MY_assert_equal ( computed , 7.467037797136421240 );
//
computed = assert_computedigits ( 3.1415926 , %pi , 2 );
-MY_assert_equal ( computed , 25.80496264389331884 );
+MY_assert_equal ( computed , 24.804962643893318841037 );
//
computed = assert_computedigits ( [0 0 1 1] , [0 1 0 1] );
MY_assert_equal ( computed , [dmax 0 0 dmax] );
MY_assert_equal ( computed , 0 );
//
computed = assert_computedigits ( 1.2345 + %i*6.7891 , 1.23456789 + %i*6.789123456 );
-MY_assert_equal ( computed , 4.259709168495138698 );
+MY_assert_equal ( computed , 3.9586791728311578886235 );
//
// The sign bit of the number of digits may be wrong because
// ieee(2); z=max(-0,0); 1/z is -%inf
computed = assert_computedigits ( 0 , 1.e-305 );
MY_assert_equal ( 1/computed , %inf );
ieee(back);
+//
+// An empirically found test case
+a = [
+3.982729777831130693D-59
+2.584939414228211484D-26
+4.391531370352049090D+43
+1.725436586898508346D+68
+];
+b = [
+3.982729777831130693D-59
+2.584939414228211484D-26
+4.391531370352048595D+43
+1.725436586898508107D+68
+];
+c = assert_computedigits ( a , b , 2 );
+e = [
+53.
+53.
+51.977632
+51.678072
+];
+assert_checkalmostequal ( c , e , 1.e-7 );
+//
+// Check that the vectorization was correct, i.e. no specific
+// case in the processing of the data is forgotten.
+//
+function pI = permInverse(p)
+ // Given the permutation p, compute the
+ // inverse permutation pI.
+ N = size(p,"*")
+ pI(p) = (1:N)'
+endfunction
+a = [
+1.234567891234567891
+1.2345678912345678
+1.23456789123456
+1.234567891234
+1.2345678912
+1.23456789
+1.234567
+1.2345
+1.23
+1.2
+1.
+0.
+%nan
+%nan
+%nan
+%inf
+%inf
+%inf
+-%inf
+-%inf
+-%inf
+0.
+0.
+-0.
+-0.
+];
+N = size(a,"*");
+for k = 1 : 10
+ mprintf("Test #%d\n",k);
+ p1 = grand(1,"prm",(1:N)');
+ p2 = grand(1,"prm",(1:N)');
+ computed = a(p1);
+ expected = a(p2);
+ d1 = assert_computedigits(computed,expected);
+ // Permute both computed and expected with the same permutation p3:
+ // d must not change.
+ p3 = grand(1,"prm",(1:N)');
+ computedP = computed(p3);
+ expectedP = expected(p3);
+ d2 = assert_computedigits(computedP,expectedP);
+ // Apply inverse permutation on d2.
+ pI = permInverse(p3);
+ d2 = d2(pI);
+ assert_checkequal(d1,d2);
+end
+Test #1
+Test #2
+Test #3
+Test #4
+Test #5
+Test #6
+Test #7
+Test #8
+Test #9
+Test #10