improve fixpoint

src/math.ml

@@ -76,7 +76,7 @@ let rec bisect_increasing ?(factor = 1e-10) ?(relax = 0.0) ?low
         if m < xfmin then false
         else
           let cf = classify_float (fx *. yfmin) in
-          (cf = FP_nan && cf = FP_infinite) || fx < yfmin
+          (cf = FP_nan || cf = FP_infinite) || fx < yfmin
       in
       let width = bmax -. bmin in
       if abs_float fx < factor || width < factor /. 10. then (*stop *) (m, fx)

src/rational_fraction.ml

@@ -224,29 +224,41 @@ module Compute (Bt : ZoneGraphInput.BoundType) (P : Polynomial.S) = struct
   let compute_z ?(verbose = 1) ?finish ?seed:_ expected_length rg =
     let w = rg.ZoneGraph.statelist.(rg.init).loc_weight in
     let dw = diff w z in
-    let ((_, denom) as size) = var z *.. dw /.. w in
+    let ((num, denom) as size) = var z *.. dw /.. w in
     if verbose > 1 then Format.printf "\nGenerating function : %a@." print size;
     let first_pole = P.find_root_sturm denom z 0.0 in
-    let tf p zv =
-      match to_float @@ apply_const p z (F.of_float zv) with
-      | Some x -> x
-      | None -> failwith "Fail to evaluate"
+    let elf = F.of_float expected_length in
+    let eq = P.( -.. ) num (P.( **. ) elf denom) in
+    let eq' = P.diff eq z in
+
+    let tf zv =
+      match
+        ( P.to_float @@ P.apply_const eq z (F.of_float zv),
+          P.to_float @@ P.apply_const eq' z (F.of_float zv) )
+      with
+      | Some x, Some y -> (x, y)
+      | _ -> failwith "Fail to evaluate"
     in
 
     let compv zval =
       if verbose > 3 then Format.printf "Size length for z: %g -> " zval;
-      let s2 = tf size zval in
-      let size = s2 -. expected_length in
+      let s2, s2' = tf zval in
 
-      if verbose > 3 then Format.printf "(%g)@." size;
-      (size, 0.0)
+      if verbose > 3 then Format.printf "(%g)@." s2;
+      (s2, s2')
     in
-    let x, v =
-      Math.bisect_increasing ~up_bound:false (0.0, first_pole) compv
+    let x =
+      Math.newton_raphson_iterate ~bound:(0.0, first_pole) compv 0.0
       (*Common.newton_raphson_iterate ~bound:(0.0, Float.infinity) ~max_iter:1000
         compv 0.1*)
     in
 
+    let v =
+      match to_float @@ apply_const size z (F.of_float x) with
+      | Some x -> x
+      | _ -> failwith "Fail to evaluate"
+    in
+
     Option.iter (fun f -> f (x, first_pole, v +. expected_length)) finish;
 
     Some x

test.t/run.t

@@ -1138,43 +1138,13 @@ Test expected size
   Computing Distribution[exact rational; RationalFraction; clocks:{}; vars:{t; z; s; }] -> 3: [] 
   Computing boltzmann parameter z for E[N]=20 
   Generating function : (-2z)/(-2+6z-4z²)
-  Size length for z: 0 -> (-20)
-  Size length for z: 0.5 -> (8.58993e+09)
-  Size length for z: 0 -> (-20)
-  Size length for z: 0.25 -> (-19.3333)
-  Size length for z: 0.375 -> (-17.6)
-  Size length for z: 0.4375 -> (-13.7778)
-  Size length for z: 0.46875 -> (-5.88235)
-  Size length for z: 0.484375 -> (10.0606)
-  Size length for z: 0.476562 -> (-0.577114)
-  Size length for z: 0.480469 -> (3.67519)
-  Size length for z: 0.478516 -> (1.35512)
-  Size length for z: 0.477539 -> (0.346851)
-  Size length for z: 0.477051 -> (-0.124998)
-  Size length for z: 0.477295 -> (0.108381)
-  Size length for z: 0.477173 -> (-0.00893477)
-  Size length for z: 0.477234 -> (0.0495651)
-  Size length for z: 0.477203 -> (0.0202759)
-  Size length for z: 0.477188 -> (0.00566074)
-  Size length for z: 0.47718 -> (-0.00163946)
-  Size length for z: 0.477184 -> (0.00201003)
-  Size length for z: 0.477182 -> (0.000185131)
-  Size length for z: 0.477181 -> (-0.000727203)
-  Size length for z: 0.477182 -> (-0.000271046)
-  Size length for z: 0.477182 -> (-4.296e-05)
-  Size length for z: 0.477182 -> (7.10847e-05)
-  Size length for z: 0.477182 -> (1.40622e-05)
-  Size length for z: 0.477182 -> (-1.4449e-05)
-  Size length for z: 0.477182 -> (-1.93392e-07)
-  Size length for z: 0.477182 -> (6.9344e-06)
-  Size length for z: 0.477182 -> (3.3705e-06)
-  Size length for z: 0.477182 -> (1.58856e-06)
-  Size length for z: 0.477182 -> (6.97582e-07)
-  Size length for z: 0.477182 -> (2.52095e-07)
-  Size length for z: 0.477182 -> (2.93518e-08)
-  Size length for z: 0.477182 -> (-8.20199e-08)
-  Size length for z: 0.477182 -> (-2.63341e-08)
-  radius of convergence:  0.5 -> z=0.477182;E[N]=20  
+  Size length for z: 0 -> (40)
+  Size length for z: 0.327869 -> (8.59984)
+  Size length for z: 0.451535 -> (1.22346)
+  Size length for z: 0.476125 -> (0.0483731)
+  Size length for z: 0.47718 -> (8.91633e-05)
+  Size length for z: 0.477182 -> (3.05182e-10)
+  radius of convergence:  0.5 -> z=0.477182;E[N]=40  
   ccaccacbccccccccbcaaaaaaaacbbbbcaa
   
   aaaccaacbbbccbbbbccbbbbbbccbbbccbccccbccccbccccccbccbcaaa