mpmath auto-tune precision

examples/bug1.cc

@@ -34,7 +34,7 @@
 int main(int argc, char *argv[]) {
   try {
     nmie::MultiLayerMieApplied<double> multi_layer_mie;
-    const std::complex<double> index_Si(4, 0.01);
+    const std::complex<double> index_Si(4, -1);
     double delta = 1e-5;
     double core_r = 100; //nm Si
     double host = 2.;

src/special-functions-impl.hpp

@@ -111,6 +111,8 @@ void evalBackwardR (const std::complex<FloatType> z,
   for (int n = nmax -1; n >= 0; n--) {
     r[n] = -static_cast<FloatType>(1)/r[n+1] + static_cast<FloatType>(2*n+1)/z;
   }
+//  r[0] = nmm::cos(z)/nmm::sin(z);
+  nmm::cout << "R0 = " << r[0] <<" at arg = "<<z<<'\n';
 }
 
 void convertRtoD1(const std::complex<FloatType> z,
@@ -138,9 +140,9 @@ void evalDownwardD1 (const std::complex<FloatType> z,
     D1[n - 1] = static_cast<FloatType>(n)*zinv - c_one/
         (D1[n] + static_cast<FloatType>(n)*zinv);
   }
-// //   r0 = cot(z)
-//  std::complex<FloatType> r0 = nmm::cos(z)/nmm::sin(z);
-//  D1[0] = r0; // - n/mx;
+ //   r0 = cot(z)
+  std::complex<FloatType> r0 = nmm::cos(z)/nmm::sin(z);
+  D1[0] = r0; // - n/mx;
 }
 
 // ********************************************************************** //

tests/special-functions-test-generator.py

@@ -17,6 +17,12 @@ Du_test = [
 [100,   [10,  -10,],  2.071124, 1.836785, 'l'],
 [10000, [10,  -10],   2.005914, 1.795393, 'm'],
 [80, [1.05,  1],   0, 0, 'Yang'],
+[1, [mp.pi,  1],   0, 0, 'pi'],
+[1, [mp.pi,  -1],   0, 0, 'pi'],
+[1, [mp.pi,  mp.pi],   0, 0, 'pi'],
+[1, [2*mp.pi,  -1],   0, 0, 'pi'],
+[1, [2*mp.pi,  mp.pi],   0, 0, 'pi'],
+[1, [2*mp.pi,  1],   0, 0, 'pi'],
 ]
 # // Dtest refractive index is m={1.05,1}, the size parameter is x = 80
 n_list = [0,1,30,50,60,70,75,80,85,90,99,116,130];
@@ -53,40 +59,72 @@ def psi(n,z):
 # to compare r(n,z) with Wolfram Alpha
 # n=49, z=1.3-2.1i,  SphericalBesselJ[n-1,z]/SphericalBesselJ[n,z]
 def r(n,z):
-    return psi(n-1,z)/psi(n,z)
+    if n > 0:
+        return psi(n-1,z)/psi(n,z)
+    return mp.cos(z)/mp.sin(z)
 
 
 def D1(n,z):
     return r(n,z) - n/z
 
 
-def main ():
-    mp.mp.dps = 30
-    output_dps = 10
-    elements_of_n_list = 15
-    out_filename = 'test_spec_functions_data.h'
-    zlist = get_z_values(Du_test)
+def get_test_data_nlist(z_record, output_dps, n):
+    x = str(z_record[0])
+    mr = str(z_record[1][0])
+    mi = str(z_record[1][1])
+    z_str = ''
+    try:
+        z = mp.mpf(x)*mp.mpc(mr, mi)
+        D1nz = D1(n,z)
+        z_str=('{{'+
+                 mp.nstr(z.real, output_dps*2)+','+
+                 mp.nstr(z.imag, output_dps*2)+'},'+
+                 str(n)+',{'+
+                 mp.nstr(D1nz.real, output_dps)+','+
+                 mp.nstr(D1nz.imag, output_dps)+'},'+
+                 mp.nstr(mp.fabs(D1nz.real* 10**-output_dps),2)+',' +
+                 mp.nstr(mp.fabs(D1nz.imag* 10**-output_dps),2)+
+                 '},\n')
+    except:
+        pass
+    return z_str
+
+def get_test_data(Du_test, output_dps, max_num_elements_of_n_list):
     output_str = ('// complex(z), n, complex(D1_n(z)), abs_err_real, abs_err_imag\n'+
-            'D1_test_'+str(output_dps)+'digits = {\n')
-    for z in zlist:
-        n_list = get_n_list(z, elements_of_n_list)
+    'std::vector< std::tuple< std::complex<double>, int, std::complex<double>, double, double > >'+
+    'D1_test_'+str(output_dps)+'digits = {\n')
+    for z_record in Du_test:
+        x = str(z_record[0])
+        mr = str(z_record[1][0])
+        mi = str(z_record[1][1])
+        mp.mp.dps = 20
+        z = mp.mpf(x)*mp.mpc(mr, mi)
+        n_list = get_n_list(z, max_num_elements_of_n_list)
         print(z, n_list)
         for n in n_list:
-            try:
-                D1nz = D1(n,z)
-            except:
-                pass
-            print(z, n)
-            output_str+=('{{'+
-                  mp.nstr(z.real, output_dps)+','+
-                  mp.nstr(z.imag, output_dps)+'},'+
-                  str(n)+',{'+
-                  mp.nstr(D1nz.real, output_dps)+','+
-                  mp.nstr(D1nz.imag, output_dps)+'},'+
-                  mp.nstr(mp.fabs(D1nz.real* 10**-output_dps),2)+',' +
-                  mp.nstr(mp.fabs(D1nz.imag* 10**-output_dps),2)+
-                  '},\n')
+            mp.mp.dps = 20
+            old_z_string = get_test_data_nlist(z_record, output_dps, n)
+            mp.mp.dps = 37
+            new_z_string = get_test_data_nlist(z_record, output_dps, n)
+            while old_z_string != new_z_string:
+                new_dps = int(mp.mp.dps * 1.41)
+                if new_dps > 300: break
+                mp.mp.dps = new_dps
+                print("New dps = ", mp.mp.dps)
+                old_z_string = new_z_string
+                new_z_string = get_test_data_nlist(z_record, output_dps, n)
+
+            output_str += new_z_string
     output_str += '};\n'
+    return output_str
+
+
+def main ():
+    output_dps = 10
+    max_num_elements_of_nlist = 15
+
+    out_filename = 'test_spec_functions_data.h'
+    output_str = get_test_data(Du_test, output_dps, max_num_elements_of_nlist)
     with open(out_filename, 'w') as out_file:
         out_file.write(output_str)
 main()

tests/test_Riccati_Bessel_logarithmic_derivative.cc

@@ -1,5 +1,6 @@
 #include "gtest/gtest.h"
 #include "../src/nmie-impl.hpp"
+#include "test_spec_functions_data.h"
 // From W. Yang APPLIED OPTICS Vol. 42, No. 9,  20 March 2003
 // Dtest refractive index is m={1.05,1}, the size parameter is x = 80
 std::vector<int> Dtest_n({0,1,30,50,60,70,75,80,85,90,99,116,130});
@@ -40,6 +41,10 @@ std::vector< std::complex<double>>
               {-0.34969,0.10437e+01 },{-0.46296,0.10809e+01 },
               {-0.56047,0.11206e+01 }});
 
+TEST(D1test, mpmath_generated_input) {
+
+}
+
 
 //TEST(D1test, DISABLED_WYang_data){
 TEST(D1test, WYang_data){

tests/test_spec_functions_data.h

@@ -1,5 +1,5 @@
 // complex(z), n, complex(D1_n(z)), abs_err_real, abs_err_imag
-D1_test_10digits = {
+std::vector< std::tuple< std::complex<double>, int, std::complex<double>, double, double > >D1_test_10digits = {
 {{0.07425,0.0},0,{13.44325437,0.0},1.3e-9,0.0},
 {{0.07425,0.0},1,{26.9211746,0.0},2.7e-9,0.0},
 {{0.07425,0.0},2,{40.39343233,0.0},4.0e-9,0.0},
@@ -63,10 +63,6 @@ D1_test_10digits = {
 {{13300.0,-0.1},167,{0.4753499708,0.1224512411},4.8e-11,1.2e-11},
 {{13300.0,-0.1},349,{0.05202548332,0.09987187722},5.2e-12,1.0e-11},
 {{13300.0,-0.1},725,{0.6178835263,0.1379545294},6.2e-11,1.4e-11},
-{{13300.0,-0.1},1508,{0.6178835263,0.1379545294},6.2e-11,1.4e-11},
-{{13300.0,-0.1},3137,{0.6178835263,0.1379545294},6.2e-11,1.4e-11},
-{{13300.0,-0.1},6522,{0.6178835263,0.1379545294},6.2e-11,1.4e-11},
-{{13300.0,-0.1},13560,{0.6178835263,0.1379545294},6.2e-11,1.4e-11},
 {{0.0825,-0.055},0,{8.364112571,5.612760191},8.4e-10,5.6e-10},
 {{0.0825,-0.055},1,{16.76671786,11.19981666},1.7e-9,1.1e-9},
 {{0.0825,-0.055},2,{25.16303989,16.7910761},2.5e-9,1.7e-9},
@@ -145,11 +141,6 @@ D1_test_10digits = {
 {{100000.0,-100000.0},238,{1.42205711e-6,1.0},1.4e-16,1.0e-10},
 {{100000.0,-100000.0},594,{8.835794179e-6,1.0},8.8e-16,1.0e-10},
 {{100000.0,-100000.0},1481,{5.48713243e-5,1.000000001},5.5e-15,1.0e-10},
-{{100000.0,-100000.0},3689,{5.48713243e-5,1.000000001},5.5e-15,1.0e-10},
-{{100000.0,-100000.0},9189,{5.48713243e-5,1.000000001},5.5e-15,1.0e-10},
-{{100000.0,-100000.0},22888,{5.48713243e-5,1.000000001},5.5e-15,1.0e-10},
-{{100000.0,-100000.0},57009,{5.48713243e-5,1.000000001},5.5e-15,1.0e-10},
-{{100000.0,-100000.0},141994,{5.48713243e-5,1.000000001},5.5e-15,1.0e-10},
 {{84.0,80.0},0,{0.0,-1.0},0.0,1.0e-10},
 {{84.0,80.0},1,{7.464603828e-5,-0.9999958865},7.5e-15,1.0e-10},
 {{84.0,80.0},2,{0.0002239401908,-0.9999876764},2.2e-14,1.0e-10},
@@ -164,4 +155,74 @@ D1_test_10digits = {
 {{84.0,80.0},77,{0.2216380379,-1.012470814},2.2e-11,1.0e-10},
 {{84.0,80.0},114,{0.4554315865,-1.075499094},4.6e-11,1.1e-10},
 {{84.0,80.0},170,{0.8526912334,-1.272209503},8.5e-11,1.3e-10},
+{{3.1415926535897932385,1.0},0,{-2.677328349e-283,-1.313035285},2.7e-293,1.3e-10},
+{{3.1415926535897932385,1.0},1,{-0.1054547331,-0.6835251095},1.1e-11,6.8e-11},
+{{3.1415926535897932385,1.0},2,{0.3756255263,-0.5129750816},3.8e-11,5.1e-11},
+{{3.1415926535897932385,1.0},3,{0.7838512831,-0.5200731151},7.8e-11,5.2e-11},
+{{3.1415926535897932385,1.0},5,{1.484043591,-0.6417333044},1.5e-10,6.4e-11},
+{{3.1415926535897932385,1.0},6,{1.808116244,-0.7189125881},1.8e-10,7.2e-11},
+{{3.1415926535897932385,1.0},7,{2.123479806,-0.8004524453},2.1e-10,8.0e-11},
+{{3.1415926535897932385,1.0},10,{3.041056951,-1.057735997},3.0e-10,1.1e-10},
+{{3.1415926535897932385,1.0},12,{3.639954628,-1.23443112},3.6e-10,1.2e-10},
+{{3.1415926535897932385,1.0},15,{4.528641422,-1.503065349},4.5e-10,1.5e-10},
+{{3.1415926535897932385,1.0},20,{5.99621535,-1.955597366},6.0e-10,2.0e-10},
+{{3.1415926535897932385,-1.0},0,{-2.677328349e-283,1.313035285},2.7e-293,1.3e-10},
+{{3.1415926535897932385,-1.0},1,{-0.1054547331,0.6835251095},1.1e-11,6.8e-11},
+{{3.1415926535897932385,-1.0},2,{0.3756255263,0.5129750816},3.8e-11,5.1e-11},
+{{3.1415926535897932385,-1.0},3,{0.7838512831,0.5200731151},7.8e-11,5.2e-11},
+{{3.1415926535897932385,-1.0},5,{1.484043591,0.6417333044},1.5e-10,6.4e-11},
+{{3.1415926535897932385,-1.0},6,{1.808116244,0.7189125881},1.8e-10,7.2e-11},
+{{3.1415926535897932385,-1.0},7,{2.123479806,0.8004524453},2.1e-10,8.0e-11},
+{{3.1415926535897932385,-1.0},10,{3.041056951,1.057735997},3.0e-10,1.1e-10},
+{{3.1415926535897932385,-1.0},12,{3.639954628,1.23443112},3.6e-10,1.2e-10},
+{{3.1415926535897932385,-1.0},15,{4.528641422,1.503065349},4.5e-10,1.5e-10},
+{{3.1415926535897932385,-1.0},20,{5.99621535,1.955597366},6.0e-10,2.0e-10},
+{{3.1415926535897932385,3.1415926535897932385},0,{-2.772406603e-285,-1.003741873},2.8e-295,1.0e-10},
+{{3.1415926535897932385,3.1415926535897932385},1,{0.05631023433,-0.984253293},5.6e-12,9.8e-11},
+{{3.1415926535897932385,3.1415926535897932385},2,{0.193283844,-0.9820447866},1.9e-11,9.8e-11},
+{{3.1415926535897932385,3.1415926535897932385},3,{0.3699265824,-1.027128175},3.7e-11,1.0e-10},
+{{3.1415926535897932385,3.1415926535897932385},4,{0.5559869777,-1.109667625},5.6e-11,1.1e-10},
+{{3.1415926535897932385,3.1415926535897932385},5,{0.741060987,-1.216012186},7.4e-11,1.2e-10},
+{{3.1415926535897932385,3.1415926535897932385},6,{0.9226942126,-1.337200281},9.2e-11,1.3e-10},
+{{3.1415926535897932385,3.1415926535897932385},8,{1.275880545,-1.605113981},1.3e-10,1.6e-10},
+{{3.1415926535897932385,3.1415926535897932385},11,{1.788056774,-2.038985764},1.8e-10,2.0e-10},
+{{3.1415926535897932385,3.1415926535897932385},14,{2.288005791,-2.49054718},2.3e-10,2.5e-10},
+{{3.1415926535897932385,3.1415926535897932385},18,{2.944407398,-3.105467936},2.9e-10,3.1e-10},
+{{3.1415926535897932385,3.1415926535897932385},23,{3.75611848,-3.88433136},3.8e-10,3.9e-10},
+{{6.283185307179586232,-1.0},0,{-1.773439591e-16,1.313035285},1.8e-26,1.3e-10},
+{{6.283185307179586232,-1.0},1,{-0.06304185944,0.7403873088},6.3e-12,7.4e-11},
+{{6.283185307179586232,-1.0},2,{0.2949375696,1.070592031},2.9e-11,1.1e-10},
+{{6.283185307179586232,-1.0},3,{-0.2986322035,0.9008011018},3.0e-11,9.0e-11},
+{{6.283185307179586232,-1.0},4,{-0.003223801583,0.4398947553},3.2e-13,4.4e-11},
+{{6.283185307179586232,-1.0},5,{0.3254847905,0.3236711082},3.3e-11,3.2e-11},
+{{6.283185307179586232,-1.0},7,{0.8216613358,0.2889872711},8.2e-11,2.9e-11},
+{{6.283185307179586232,-1.0},9,{1.227344826,0.3098261586},1.2e-10,3.1e-11},
+{{6.283185307179586232,-1.0},12,{1.773455192,0.3647272463},1.8e-10,3.6e-11},
+{{6.283185307179586232,-1.0},16,{2.453967796,0.4514019205},2.5e-10,4.5e-11},
+{{6.283185307179586232,-1.0},20,{3.110719695,0.543555894},3.1e-10,5.4e-11},
+{{6.283185307179586232,-1.0},27,{4.234797254,0.7099061556},4.2e-10,7.1e-11},
+{{6.283185307179586232,3.1415926535897932385},0,{-1.83641862e-18,-1.003741873},1.8e-28,1.0e-10},
+{{6.283185307179586232,3.1415926535897932385},1,{0.01415314238,-0.9809158089},1.4e-12,9.8e-11},
+{{6.283185307179586232,3.1415926535897932385},2,{0.05114736247,-0.9580400081},5.1e-12,9.6e-11},
+{{6.283185307179586232,3.1415926535897932385},3,{0.09211316558,-0.9082339054},9.2e-12,9.1e-11},
+{{6.283185307179586232,3.1415926535897932385},4,{0.1846022792,-0.8324588456},1.8e-11,8.3e-11},
+{{6.283185307179586232,3.1415926535897932385},6,{0.4906578918,-0.7594698541},4.9e-11,7.6e-11},
+{{6.283185307179586232,3.1415926535897932385},7,{0.6556002657,-0.7661759617},6.6e-11,7.7e-11},
+{{6.283185307179586232,3.1415926535897932385},10,{1.126363789,-0.8654702596},1.1e-10,8.7e-11},
+{{6.283185307179586232,3.1415926535897932385},13,{1.564605169,-1.013646818},1.6e-10,1.0e-10},
+{{6.283185307179586232,3.1415926535897932385},17,{2.121175275,-1.23751537},2.1e-10,1.2e-10},
+{{6.283185307179586232,3.1415926535897932385},22,{2.794289612,-1.534314595},2.8e-10,1.5e-10},
+{{6.283185307179586232,3.1415926535897932385},28,{3.585638859,-1.901082509},3.6e-10,1.9e-10},
+{{6.283185307179586232,1.0},0,{-1.773439591e-16,-1.313035285},1.8e-26,1.3e-10},
+{{6.283185307179586232,1.0},1,{-0.06304185944,-0.7403873088},6.3e-12,7.4e-11},
+{{6.283185307179586232,1.0},2,{0.2949375696,-1.070592031},2.9e-11,1.1e-10},
+{{6.283185307179586232,1.0},3,{-0.2986322035,-0.9008011018},3.0e-11,9.0e-11},
+{{6.283185307179586232,1.0},4,{-0.003223801583,-0.4398947553},3.2e-13,4.4e-11},
+{{6.283185307179586232,1.0},5,{0.3254847905,-0.3236711082},3.3e-11,3.2e-11},
+{{6.283185307179586232,1.0},7,{0.8216613358,-0.2889872711},8.2e-11,2.9e-11},
+{{6.283185307179586232,1.0},9,{1.227344826,-0.3098261586},1.2e-10,3.1e-11},
+{{6.283185307179586232,1.0},12,{1.773455192,-0.3647272463},1.8e-10,3.6e-11},
+{{6.283185307179586232,1.0},16,{2.453967796,-0.4514019205},2.5e-10,4.5e-11},
+{{6.283185307179586232,1.0},20,{3.110719695,-0.543555894},3.1e-10,5.4e-11},
+{{6.283185307179586232,1.0},27,{4.234797254,-0.7099061556},4.2e-10,7.1e-11},
 };