#BEGIN OUTFILE1 # Begin Function number 3 display_alot := proc(iter) global glob_iolevel, INFO, DEBUGMASSIVE, glob_max_terms, ALWAYS, DEBUGL, #Top Generate Globals Decl glob_current_iter, glob_hmin_init, years_in_century, hours_in_day, glob_display_flag, glob_max_trunc_err, glob_initial_pass, glob_almost_1, glob_curr_iter_when_opt, glob_max_sec, glob_optimal_start, glob_dump_analytic, glob_not_yet_finished, djd_debug, glob_optimal_expect_sec, glob_warned, glob_optimal_clock_start_sec, glob_reached_optimal_h, centuries_in_millinium, glob_warned2, glob_disp_incr, glob_clock_sec, glob_log10normmin, glob_iter, glob_no_eqs, min_in_hour, glob_subiter_method, glob_normmax, glob_start, glob_orig_start_sec, glob_max_hours, glob_large_float, glob_hmin, days_in_year, djd_debug2, glob_smallish_float, glob_small_float, glob_max_iter, glob_look_poles, glob_h, glob_dump, glob_max_minutes, glob_log10abserr, glob_abserr, glob_log10_abserr, glob_optimal_done, glob_clock_start_sec, glob_max_opt_iter, glob_percent_done, glob_unchanged_h_cnt, glob_max_rel_trunc_err, glob_hmax, glob_not_yet_start_msg, sec_in_min, glob_log10relerr, MAX_UNCHANGED, glob_relerr, glob_log10_relerr, glob_last_good_h, glob_html_log, #Bottom Generate Globals Decl #BEGIN CONST array_const_1D0, array_const_1, array_const_0D0, #END CONST array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_y_init, array_1st_rel_error, array_y, array_x, array_m1, array_last_rel_error, array_pole, array_norms, array_type_pole, array_y_higher_work2, array_real_pole, array_y_higher, array_y_set_initial, array_y_higher_work, array_complex_pole, array_poles, glob_last; local abserr, analytic_val_y, ind_var, numeric_val, relerr, term_no; #TOP DISPLAY ALOT if (iter >= 0) then # if number 1 ind_var := array_x[1]; omniout_float(ALWAYS,"x[1] ",33,ind_var,20," "); analytic_val_y := exact_soln_y(ind_var); omniout_float(ALWAYS,"y[1] (analytic) ",33,analytic_val_y,20," "); term_no := 1; numeric_val := array_y[term_no]; abserr := abs(numeric_val - analytic_val_y); omniout_float(ALWAYS,"y[1] (numeric) ",33,numeric_val,20," "); if (abs(analytic_val_y) <> 0.0) then # if number 2 relerr := abserr*100.0/abs(analytic_val_y); else relerr := -1.0 ; fi;# end if 2 ; if glob_iter = 1 then # if number 2 array_1st_rel_error[1] := relerr; else array_last_rel_error[1] := relerr; fi;# end if 2 ; omniout_float(ALWAYS,"absolute error ",4,abserr,20," "); omniout_float(ALWAYS,"relative error ",4,relerr,20,"%"); omniout_float(ALWAYS,"h ",4,glob_h,20," "); #BOTTOM DISPLAY ALOT fi;# end if 1 ; # End Function number 3 end; # Begin Function number 4 adjust_for_pole := proc(h_param) global glob_iolevel, INFO, DEBUGMASSIVE, glob_max_terms, ALWAYS, DEBUGL, #Top Generate Globals Decl glob_current_iter, glob_hmin_init, years_in_century, hours_in_day, glob_display_flag, glob_max_trunc_err, glob_initial_pass, glob_almost_1, glob_curr_iter_when_opt, glob_max_sec, glob_optimal_start, glob_dump_analytic, glob_not_yet_finished, djd_debug, glob_optimal_expect_sec, glob_warned, glob_optimal_clock_start_sec, glob_reached_optimal_h, centuries_in_millinium, glob_warned2, glob_disp_incr, glob_clock_sec, glob_log10normmin, glob_iter, glob_no_eqs, min_in_hour, glob_subiter_method, glob_normmax, glob_start, glob_orig_start_sec, glob_max_hours, glob_large_float, glob_hmin, days_in_year, djd_debug2, glob_smallish_float, glob_small_float, glob_max_iter, glob_look_poles, glob_h, glob_dump, glob_max_minutes, glob_log10abserr, glob_abserr, glob_log10_abserr, glob_optimal_done, glob_clock_start_sec, glob_max_opt_iter, glob_percent_done, glob_unchanged_h_cnt, glob_max_rel_trunc_err, glob_hmax, glob_not_yet_start_msg, sec_in_min, glob_log10relerr, MAX_UNCHANGED, glob_relerr, glob_log10_relerr, glob_last_good_h, glob_html_log, #Bottom Generate Globals Decl #BEGIN CONST array_const_1D0, array_const_1, array_const_0D0, #END CONST array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_y_init, array_1st_rel_error, array_y, array_x, array_m1, array_last_rel_error, array_pole, array_norms, array_type_pole, array_y_higher_work2, array_real_pole, array_y_higher, array_y_set_initial, array_y_higher_work, array_complex_pole, array_poles, glob_last; local hnew, sz2, tmp; #TOP ADJUST FOR POLE hnew := h_param; glob_normmax := glob_small_float; if (abs(array_y_higher[1,1]) > glob_small_float) then # if number 1 tmp := abs(array_y_higher[1,1]); if (tmp < glob_normmax) then # if number 2 glob_normmax := tmp; fi;# end if 2 fi;# end if 1 ; if (glob_look_poles and (abs(array_pole[1]) > glob_small_float) and (array_pole[1] <> glob_large_float)) then # if number 1 sz2 := array_pole[1]/10.0; if (sz2 < hnew) then # if number 2 omniout_float(INFO,"glob_h adjusted to ",20,h_param,12,"due to singularity."); omniout_str(INFO,"Reached Optimal"); newline(); return(hnew); fi;# end if 2 fi;# end if 1 ; if (not glob_reached_optimal_h) then # if number 1 glob_reached_optimal_h := true; glob_curr_iter_when_opt := glob_current_iter; glob_optimal_clock_start_sec := elapsed_time_seconds(); glob_optimal_start := array_x[1]; fi;# end if 1 ; hnew := sz2; #END block #BOTTOM ADJUST FOR POLE # End Function number 4 end; # Begin Function number 5 prog_report := proc(x_start,x_end) global glob_iolevel, INFO, DEBUGMASSIVE, glob_max_terms, ALWAYS, DEBUGL, #Top Generate Globals Decl glob_current_iter, glob_hmin_init, years_in_century, hours_in_day, glob_display_flag, glob_max_trunc_err, glob_initial_pass, glob_almost_1, glob_curr_iter_when_opt, glob_max_sec, glob_optimal_start, glob_dump_analytic, glob_not_yet_finished, djd_debug, glob_optimal_expect_sec, glob_warned, glob_optimal_clock_start_sec, glob_reached_optimal_h, centuries_in_millinium, glob_warned2, glob_disp_incr, glob_clock_sec, glob_log10normmin, glob_iter, glob_no_eqs, min_in_hour, glob_subiter_method, glob_normmax, glob_start, glob_orig_start_sec, glob_max_hours, glob_large_float, glob_hmin, days_in_year, djd_debug2, glob_smallish_float, glob_small_float, glob_max_iter, glob_look_poles, glob_h, glob_dump, glob_max_minutes, glob_log10abserr, glob_abserr, glob_log10_abserr, glob_optimal_done, glob_clock_start_sec, glob_max_opt_iter, glob_percent_done, glob_unchanged_h_cnt, glob_max_rel_trunc_err, glob_hmax, glob_not_yet_start_msg, sec_in_min, glob_log10relerr, MAX_UNCHANGED, glob_relerr, glob_log10_relerr, glob_last_good_h, glob_html_log, #Bottom Generate Globals Decl #BEGIN CONST array_const_1D0, array_const_1, array_const_0D0, #END CONST array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_y_init, array_1st_rel_error, array_y, array_x, array_m1, array_last_rel_error, array_pole, array_norms, array_type_pole, array_y_higher_work2, array_real_pole, array_y_higher, array_y_set_initial, array_y_higher_work, array_complex_pole, array_poles, glob_last; local clock_sec, opt_clock_sec, clock_sec1, expect_sec, left_sec, percent_done, total_clock_sec; #TOP PROGRESS REPORT clock_sec1 := elapsed_time_seconds(); total_clock_sec := convfloat(clock_sec1) - convfloat(glob_orig_start_sec); glob_clock_sec := convfloat(clock_sec1) - convfloat(glob_clock_start_sec); left_sec := convfloat(glob_max_sec) + convfloat(glob_orig_start_sec) - convfloat(clock_sec1); expect_sec := comp_expect_sec(convfloat(x_end),convfloat(x_start),convfloat(array_x[1]) + convfloat(glob_h) ,convfloat( clock_sec1) - convfloat(glob_orig_start_sec)); opt_clock_sec := convfloat( clock_sec1) - convfloat(glob_optimal_clock_start_sec); glob_optimal_expect_sec := comp_expect_sec(convfloat(x_end),convfloat(x_start),convfloat(array_x[1]) +convfloat( glob_h) ,convfloat( opt_clock_sec)); percent_done := comp_percent(convfloat(x_end),convfloat(x_start),convfloat(array_x[1]) + convfloat(glob_h)); glob_percent_done := percent_done; omniout_str_noeol(INFO,"Total Elapsed Time "); omniout_timestr(convfloat(total_clock_sec)); omniout_str_noeol(INFO,"Elapsed Time(since restart) "); omniout_timestr(convfloat(glob_clock_sec)); if convfloat(percent_done) < convfloat(100.0) then # if number 1 omniout_str_noeol(INFO,"Expected Time Remaining "); omniout_timestr(convfloat(expect_sec)); omniout_str_noeol(INFO,"Optimized Time Remaining "); omniout_timestr(convfloat(glob_optimal_expect_sec)); fi;# end if 1 ; omniout_str_noeol(INFO,"Time to Timeout "); omniout_timestr(convfloat(left_sec)); omniout_float(INFO, "Percent Done ",33,percent_done,4,"%"); #BOTTOM PROGRESS REPORT # End Function number 5 end; # Begin Function number 6 check_for_pole := proc() global glob_iolevel, INFO, DEBUGMASSIVE, glob_max_terms, ALWAYS, DEBUGL, #Top Generate Globals Decl glob_current_iter, glob_hmin_init, years_in_century, hours_in_day, glob_display_flag, glob_max_trunc_err, glob_initial_pass, glob_almost_1, glob_curr_iter_when_opt, glob_max_sec, glob_optimal_start, glob_dump_analytic, glob_not_yet_finished, djd_debug, glob_optimal_expect_sec, glob_warned, glob_optimal_clock_start_sec, glob_reached_optimal_h, centuries_in_millinium, glob_warned2, glob_disp_incr, glob_clock_sec, glob_log10normmin, glob_iter, glob_no_eqs, min_in_hour, glob_subiter_method, glob_normmax, glob_start, glob_orig_start_sec, glob_max_hours, glob_large_float, glob_hmin, days_in_year, djd_debug2, glob_smallish_float, glob_small_float, glob_max_iter, glob_look_poles, glob_h, glob_dump, glob_max_minutes, glob_log10abserr, glob_abserr, glob_log10_abserr, glob_optimal_done, glob_clock_start_sec, glob_max_opt_iter, glob_percent_done, glob_unchanged_h_cnt, glob_max_rel_trunc_err, glob_hmax, glob_not_yet_start_msg, sec_in_min, glob_log10relerr, MAX_UNCHANGED, glob_relerr, glob_log10_relerr, glob_last_good_h, glob_html_log, #Bottom Generate Globals Decl #BEGIN CONST array_const_1D0, array_const_1, array_const_0D0, #END CONST array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_y_init, array_1st_rel_error, array_y, array_x, array_m1, array_last_rel_error, array_pole, array_norms, array_type_pole, array_y_higher_work2, array_real_pole, array_y_higher, array_y_set_initial, array_y_higher_work, array_complex_pole, array_poles, glob_last; local cnt, dr1, dr2, ds1, ds2, hdrc, m, n, nr1, nr2, ord_no, rad_c, rcs, rm0, rm1, rm2, rm3, rm4, found; #TOP CHECK FOR POLE #IN RADII REAL EQ = 1 #Computes radius of convergence and r_order of pole from 3 adjacent Taylor series terms. EQUATUON NUMBER 1 #Applies to pole of arbitrary r_order on the real axis, #Due to Prof. George Corliss. n := glob_max_terms; m := n - 1 - 1; while ((m >= 10) and ((abs(array_y_higher[1,m]) < glob_small_float) or (abs(array_y_higher[1,m-1]) < glob_small_float) or (abs(array_y_higher[1,m-2]) < glob_small_float ))) do # do number 2 m := m - 1; od;# end do number 2 ; if (m > 10) then # if number 1 rm0 := array_y_higher[1,m]/array_y_higher[1,m-1]; rm1 := array_y_higher[1,m-1]/array_y_higher[1,m-2]; hdrc := convfloat(m-1)*rm0-convfloat(m-2)*rm1; if (abs(hdrc) > glob_small_float) then # if number 2 rcs := glob_h/hdrc; ord_no := convfloat(m-1)*rm0/hdrc - convfloat(m) + 2.0; array_real_pole[1,1] := rcs; array_real_pole[1,2] := ord_no; else array_real_pole[1,1] := glob_large_float; array_real_pole[1,2] := glob_large_float; fi;# end if 2 else array_real_pole[1,1] := glob_large_float; array_real_pole[1,2] := glob_large_float; fi;# end if 1 ; #BOTTOM RADII REAL EQ = 1 #TOP RADII COMPLEX EQ = 1 #Computes radius of convergence for complex conjugate pair of poles. #from 6 adjacent Taylor series terms #Also computes r_order of poles. #Due to Manuel Prieto. #With a correction by Dennis J. Darland n := glob_max_terms - 1 - 1; cnt := 0; while ((cnt < 5) and (n >= 10)) do # do number 2 if (abs(array_y_higher[1,n]) > glob_small_float) then # if number 1 cnt := cnt + 1; else cnt := 0; fi;# end if 1 ; n := n - 1; od;# end do number 2 ; m := n + cnt; if (m <= 10) then # if number 1 array_complex_pole[1,1] := glob_large_float; array_complex_pole[1,2] := glob_large_float; elif (abs(array_y_higher[1,m]) >= (glob_large_float)) or (abs(array_y_higher[1,m-1]) >=(glob_large_float)) or (abs(array_y_higher[1,m-2]) >= (glob_large_float)) or (abs(array_y_higher[1,m-3]) >= (glob_large_float)) or (abs(array_y_higher[1,m-4]) >= (glob_large_float)) or (abs(array_y_higher[1,m-5]) >= (glob_large_float)) then # if number 2 array_complex_pole[1,1] := glob_large_float; array_complex_pole[1,2] := glob_large_float; else rm0 := (array_y_higher[1,m])/(array_y_higher[1,m-1]); rm1 := (array_y_higher[1,m-1])/(array_y_higher[1,m-2]); rm2 := (array_y_higher[1,m-2])/(array_y_higher[1,m-3]); rm3 := (array_y_higher[1,m-3])/(array_y_higher[1,m-4]); rm4 := (array_y_higher[1,m-4])/(array_y_higher[1,m-5]); nr1 := convfloat(m-1)*rm0 - 2.0*convfloat(m-2)*rm1 + convfloat(m-3)*rm2; nr2 := convfloat(m-2)*rm1 - 2.0*convfloat(m-3)*rm2 + convfloat(m-4)*rm3; dr1 := (-1.0)/rm1 + 2.0/rm2 - 1.0/rm3; dr2 := (-1.0)/rm2 + 2.0/rm3 - 1.0/rm4; ds1 := 3.0/rm1 - 8.0/rm2 + 5.0/rm3; ds2 := 3.0/rm2 - 8.0/rm3 + 5.0/rm4; if ((abs(nr1 * dr2 - nr2 * dr1) <= glob_small_float) or (abs(dr1) <= glob_small_float)) then # if number 3 array_complex_pole[1,1] := glob_large_float; array_complex_pole[1,2] := glob_large_float; else if (abs(nr1*dr2 - nr2 * dr1) > glob_small_float) then # if number 4 rcs := ((ds1*dr2 - ds2*dr1 +dr1*dr2)/(nr1*dr2 - nr2 * dr1)); #(Manuels) rcs := (ds1*dr2 - ds2*dr1)/(nr1*dr2 - nr2 * dr1) ord_no := (rcs*nr1 - ds1)/(2.0*dr1) -convfloat(m)/2.0; if (abs(rcs) > glob_small_float) then # if number 5 if (rcs > 0.0) then # if number 6 rad_c := sqrt(rcs) * glob_h; else rad_c := glob_large_float; fi;# end if 6 else rad_c := glob_large_float; ord_no := glob_large_float; fi;# end if 5 else rad_c := glob_large_float; ord_no := glob_large_float; fi;# end if 4 fi;# end if 3 ; array_complex_pole[1,1] := rad_c; array_complex_pole[1,2] := ord_no; fi;# end if 2 ; #BOTTOM RADII COMPLEX EQ = 1 found := false; #TOP WHICH RADII EQ = 1 if not found and ((array_real_pole[1,1] = glob_large_float) or (array_real_pole[1,2] = glob_large_float)) and ((array_complex_pole[1,1] <> glob_large_float) and (array_complex_pole[1,2] <> glob_large_float)) and ((array_complex_pole[1,1] > 0.0) and (array_complex_pole[1,2] > 0.0)) then # if number 2 array_poles[1,1] := array_complex_pole[1,1]; array_poles[1,2] := array_complex_pole[1,2]; found := true; array_type_pole[1] := 2; if (glob_display_flag) then # if number 3 omniout_str(ALWAYS,"Complex estimate of poles used"); fi;# end if 3 ; fi;# end if 2 ; if not found and ((array_real_pole[1,1] <> glob_large_float) and (array_real_pole[1,2] <> glob_large_float) and (array_real_pole[1,1] > 0.0) and (array_real_pole[1,2] > 0.0) and ((array_complex_pole[1,1] = glob_large_float) or (array_complex_pole[1,2] = glob_large_float) or (array_complex_pole[1,1] <= 0.0 ) or (array_complex_pole[1,2] <= 0.0))) then # if number 2 array_poles[1,1] := array_real_pole[1,1]; array_poles[1,2] := array_real_pole[1,2]; found := true; array_type_pole[1] := 1; if (glob_display_flag) then # if number 3 omniout_str(ALWAYS,"Real estimate of pole used"); fi;# end if 3 ; fi;# end if 2 ; if not found and (((array_real_pole[1,1] = glob_large_float) or (array_real_pole[1,2] = glob_large_float)) and ((array_complex_pole[1,1] = glob_large_float) or (array_complex_pole[1,2] = glob_large_float))) then # if number 2 array_poles[1,1] := glob_large_float; array_poles[1,2] := glob_large_float; found := true; array_type_pole[1] := 3; if (glob_display_flag) then # if number 3 omniout_str(ALWAYS,"NO POLE"); fi;# end if 3 ; fi;# end if 2 ; if not found and ((array_real_pole[1,1] < array_complex_pole[1,1]) and (array_real_pole[1,1] > 0.0) and (array_real_pole[1,2] > 0.0)) then # if number 2 array_poles[1,1] := array_real_pole[1,1]; array_poles[1,2] := array_real_pole[1,2]; found := true; array_type_pole[1] := 1; if (glob_display_flag) then # if number 3 omniout_str(ALWAYS,"Real estimate of pole used"); fi;# end if 3 ; fi;# end if 2 ; if not found and ((array_complex_pole[1,1] <> glob_large_float) and (array_complex_pole[1,2] <> glob_large_float) and (array_complex_pole[1,1] > 0.0) and (array_complex_pole[1,2] > 0.0)) then # if number 2 array_poles[1,1] := array_complex_pole[1,1]; array_poles[1,2] := array_complex_pole[1,2]; array_type_pole[1] := 2; found := true; if (glob_display_flag) then # if number 3 omniout_str(ALWAYS,"Complex estimate of poles used"); fi;# end if 3 ; fi;# end if 2 ; if not found then # if number 2 array_poles[1,1] := glob_large_float; array_poles[1,2] := glob_large_float; array_type_pole[1] := 3; if (glob_display_flag) then # if number 3 omniout_str(ALWAYS,"NO POLE"); fi;# end if 3 ; fi;# end if 2 ; #BOTTOM WHICH RADII EQ = 1 array_pole[1] := glob_large_float; array_pole[2] := glob_large_float; #TOP WHICH RADIUS EQ = 1 if array_pole[1] > array_poles[1,1] then # if number 2 array_pole[1] := array_poles[1,1]; array_pole[2] := array_poles[1,2]; fi;# end if 2 ; #BOTTOM WHICH RADIUS EQ = 1 #BOTTOM CHECK FOR POLE display_pole(); # End Function number 6 end; # Begin Function number 7 get_norms := proc() global glob_iolevel, INFO, DEBUGMASSIVE, glob_max_terms, ALWAYS, DEBUGL, #Top Generate Globals Decl glob_current_iter, glob_hmin_init, years_in_century, hours_in_day, glob_display_flag, glob_max_trunc_err, glob_initial_pass, glob_almost_1, glob_curr_iter_when_opt, glob_max_sec, glob_optimal_start, glob_dump_analytic, glob_not_yet_finished, djd_debug, glob_optimal_expect_sec, glob_warned, glob_optimal_clock_start_sec, glob_reached_optimal_h, centuries_in_millinium, glob_warned2, glob_disp_incr, glob_clock_sec, glob_log10normmin, glob_iter, glob_no_eqs, min_in_hour, glob_subiter_method, glob_normmax, glob_start, glob_orig_start_sec, glob_max_hours, glob_large_float, glob_hmin, days_in_year, djd_debug2, glob_smallish_float, glob_small_float, glob_max_iter, glob_look_poles, glob_h, glob_dump, glob_max_minutes, glob_log10abserr, glob_abserr, glob_log10_abserr, glob_optimal_done, glob_clock_start_sec, glob_max_opt_iter, glob_percent_done, glob_unchanged_h_cnt, glob_max_rel_trunc_err, glob_hmax, glob_not_yet_start_msg, sec_in_min, glob_log10relerr, MAX_UNCHANGED, glob_relerr, glob_log10_relerr, glob_last_good_h, glob_html_log, #Bottom Generate Globals Decl #BEGIN CONST array_const_1D0, array_const_1, array_const_0D0, #END CONST array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_y_init, array_1st_rel_error, array_y, array_x, array_m1, array_last_rel_error, array_pole, array_norms, array_type_pole, array_y_higher_work2, array_real_pole, array_y_higher, array_y_set_initial, array_y_higher_work, array_complex_pole, array_poles, glob_last; local iii; if (not glob_initial_pass) then # if number 2 set_z(array_norms,glob_max_terms+1); #TOP GET NORMS iii := 1; while (iii <= glob_max_terms) do # do number 2 if (abs(array_y[iii]) > array_norms[iii]) then # if number 3 array_norms[iii] := abs(array_y[iii]); fi;# end if 3 ; iii := iii + 1; od;# end do number 2 #GET NORMS ; fi;# end if 2 ; # End Function number 7 end; # Begin Function number 8 atomall := proc() global glob_iolevel, INFO, DEBUGMASSIVE, glob_max_terms, ALWAYS, DEBUGL, #Top Generate Globals Decl glob_current_iter, glob_hmin_init, years_in_century, hours_in_day, glob_display_flag, glob_max_trunc_err, glob_initial_pass, glob_almost_1, glob_curr_iter_when_opt, glob_max_sec, glob_optimal_start, glob_dump_analytic, glob_not_yet_finished, djd_debug, glob_optimal_expect_sec, glob_warned, glob_optimal_clock_start_sec, glob_reached_optimal_h, centuries_in_millinium, glob_warned2, glob_disp_incr, glob_clock_sec, glob_log10normmin, glob_iter, glob_no_eqs, min_in_hour, glob_subiter_method, glob_normmax, glob_start, glob_orig_start_sec, glob_max_hours, glob_large_float, glob_hmin, days_in_year, djd_debug2, glob_smallish_float, glob_small_float, glob_max_iter, glob_look_poles, glob_h, glob_dump, glob_max_minutes, glob_log10abserr, glob_abserr, glob_log10_abserr, glob_optimal_done, glob_clock_start_sec, glob_max_opt_iter, glob_percent_done, glob_unchanged_h_cnt, glob_max_rel_trunc_err, glob_hmax, glob_not_yet_start_msg, sec_in_min, glob_log10relerr, MAX_UNCHANGED, glob_relerr, glob_log10_relerr, glob_last_good_h, glob_html_log, #Bottom Generate Globals Decl #BEGIN CONST array_const_1D0, array_const_1, array_const_0D0, #END CONST array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_y_init, array_1st_rel_error, array_y, array_x, array_m1, array_last_rel_error, array_pole, array_norms, array_type_pole, array_y_higher_work2, array_real_pole, array_y_higher, array_y_set_initial, array_y_higher_work, array_complex_pole, array_poles, glob_last; local kkk, order_d, adj2, temporary, term; #TOP ATOMALL #END OUTFILE1 #BEGIN ATOMHDR1 # emit pre mult $eq_no = 1 i = 1 array_tmp1[1] := (array_x[1] * (array_x[1])); #emit pre add $eq_no = 1 i = 1 array_tmp2[1] := array_tmp1[1] + array_const_1D0[1]; #emit pre div $eq_no = 1 i = 1 array_tmp3[1] := (array_const_1D0[1] / (array_tmp2[1])); #emit pre add $eq_no = 1 i = 1 array_tmp4[1] := array_const_0D0[1] + array_tmp3[1]; #emit pre assign xxx $eq_no = 1 i = 1 $min_hdrs = 5 if not array_y_set_initial[1,2] then # if number 1 if (1 <= glob_max_terms) then # if number 2 temporary := array_tmp4[1] * (glob_h ^ (1)) * factorial_3(0,1); array_y[2] := temporary; array_y_higher[1,2] := temporary; temporary := temporary / glob_h * (2.0); array_y_higher[2,1] := temporary ; fi;# end if 2 ; fi;# end if 1 ; kkk := 2; #END ATOMHDR1 #BEGIN ATOMHDR2 # emit pre mult $eq_no = 1 i = 2 array_tmp1[2] := ats(2,array_x,array_x,1); #emit pre add $eq_no = 1 i = 2 array_tmp2[2] := array_tmp1[2] + array_const_1D0[2]; #emit pre div $eq_no = 1 i = 2 array_tmp3[2] := ((array_const_1D0[2] - ats(2,array_tmp2,array_tmp3,2))/array_tmp2[1]); #emit pre add $eq_no = 1 i = 2 array_tmp4[2] := array_const_0D0[2] + array_tmp3[2]; #emit pre assign xxx $eq_no = 1 i = 2 $min_hdrs = 5 if not array_y_set_initial[1,3] then # if number 1 if (2 <= glob_max_terms) then # if number 2 temporary := array_tmp4[2] * (glob_h ^ (1)) * factorial_3(1,2); array_y[3] := temporary; array_y_higher[1,3] := temporary; temporary := temporary / glob_h * (2.0); array_y_higher[2,2] := temporary ; fi;# end if 2 ; fi;# end if 1 ; kkk := 3; #END ATOMHDR2 #BEGIN ATOMHDR3 # emit pre mult $eq_no = 1 i = 3 array_tmp1[3] := ats(3,array_x,array_x,1); #emit pre add $eq_no = 1 i = 3 array_tmp2[3] := array_tmp1[3] + array_const_1D0[3]; #emit pre div $eq_no = 1 i = 3 array_tmp3[3] := ((array_const_1D0[3] - ats(3,array_tmp2,array_tmp3,2))/array_tmp2[1]); #emit pre add $eq_no = 1 i = 3 array_tmp4[3] := array_const_0D0[3] + array_tmp3[3]; #emit pre assign xxx $eq_no = 1 i = 3 $min_hdrs = 5 if not array_y_set_initial[1,4] then # if number 1 if (3 <= glob_max_terms) then # if number 2 temporary := array_tmp4[3] * (glob_h ^ (1)) * factorial_3(2,3); array_y[4] := temporary; array_y_higher[1,4] := temporary; temporary := temporary / glob_h * (2.0); array_y_higher[2,3] := temporary ; fi;# end if 2 ; fi;# end if 1 ; kkk := 4; #END ATOMHDR3 #BEGIN ATOMHDR4 # emit pre mult $eq_no = 1 i = 4 array_tmp1[4] := ats(4,array_x,array_x,1); #emit pre add $eq_no = 1 i = 4 array_tmp2[4] := array_tmp1[4] + array_const_1D0[4]; #emit pre div $eq_no = 1 i = 4 array_tmp3[4] := ((array_const_1D0[4] - ats(4,array_tmp2,array_tmp3,2))/array_tmp2[1]); #emit pre add $eq_no = 1 i = 4 array_tmp4[4] := array_const_0D0[4] + array_tmp3[4]; #emit pre assign xxx $eq_no = 1 i = 4 $min_hdrs = 5 if not array_y_set_initial[1,5] then # if number 1 if (4 <= glob_max_terms) then # if number 2 temporary := array_tmp4[4] * (glob_h ^ (1)) * factorial_3(3,4); array_y[5] := temporary; array_y_higher[1,5] := temporary; temporary := temporary / glob_h * (2.0); array_y_higher[2,4] := temporary ; fi;# end if 2 ; fi;# end if 1 ; kkk := 5; #END ATOMHDR4 #BEGIN ATOMHDR5 # emit pre mult $eq_no = 1 i = 5 array_tmp1[5] := ats(5,array_x,array_x,1); #emit pre add $eq_no = 1 i = 5 array_tmp2[5] := array_tmp1[5] + array_const_1D0[5]; #emit pre div $eq_no = 1 i = 5 array_tmp3[5] := ((array_const_1D0[5] - ats(5,array_tmp2,array_tmp3,2))/array_tmp2[1]); #emit pre add $eq_no = 1 i = 5 array_tmp4[5] := array_const_0D0[5] + array_tmp3[5]; #emit pre assign xxx $eq_no = 1 i = 5 $min_hdrs = 5 if not array_y_set_initial[1,6] then # if number 1 if (5 <= glob_max_terms) then # if number 2 temporary := array_tmp4[5] * (glob_h ^ (1)) * factorial_3(4,5); array_y[6] := temporary; array_y_higher[1,6] := temporary; temporary := temporary / glob_h * (2.0); array_y_higher[2,5] := temporary ; fi;# end if 2 ; fi;# end if 1 ; kkk := 6; #END ATOMHDR5 #BEGIN OUTFILE3 #Top Atomall While Loop-- outfile3 while (kkk <= glob_max_terms) do # do number 1 #END OUTFILE3 #BEGIN OUTFILE4 #emit mult $eq_no = 1 array_tmp1[kkk] := ats(kkk,array_x,array_x,1); #emit add $eq_no = 1 array_tmp2[kkk] := array_tmp1[kkk] + array_const_1D0[kkk]; #emit div $eq_no = 1 array_tmp3[kkk] := ((array_const_1D0[kkk] - ats(kkk,array_tmp2,array_tmp3,2))/array_tmp2[1]); #emit add $eq_no = 1 array_tmp4[kkk] := array_const_0D0[kkk] + array_tmp3[kkk]; #emit assign $eq_no = 1 order_d := 1; if (kkk + order_d + 1 <= glob_max_terms) then # if number 1 if not array_y_set_initial[1,kkk + order_d] then # if number 2 temporary := array_tmp4[kkk] * (glob_h ^ (order_d)) / factorial_3((kkk - 1),(kkk + order_d - 1)); array_y[kkk + order_d] := temporary; array_y_higher[1,kkk + order_d] := temporary; term := kkk + order_d - 1; adj2 := 2; while (adj2 <= order_d + 1) and (term >= 1) do # do number 2 temporary := temporary / glob_h * convfp(adj2); array_y_higher[adj2,term] := temporary; adj2 := adj2 + 1; term := term - 1; od;# end do number 2 fi;# end if 2 fi;# end if 1 ; kkk := kkk + 1; od;# end do number 1 ; #BOTTOM ATOMALL #END OUTFILE4 #BEGIN OUTFILE5 # End Function number 8 end; #BEGIN ATS LIBRARY BLOCK omniout_str := proc(iolevel,str) global glob_iolevel; if (glob_iolevel >= iolevel) then printf("%s\n",str); fi; # End Function number 1 end; omniout_str_noeol := proc(iolevel,str) global glob_iolevel; if (glob_iolevel >= iolevel) then printf("%s",str); fi; # End Function number 1 end; omniout_labstr := proc(iolevel,label,str) global glob_iolevel; if (glob_iolevel >= iolevel) then print(label,str); fi; # End Function number 1 end; omniout_float := proc(iolevel,prelabel,prelen,value,vallen,postlabel) global glob_iolevel; if (glob_iolevel >= iolevel) then if vallen = 4 then printf("%-30s = %-42.4g %s \n",prelabel,value, postlabel); else printf("%-30s = %-42.32g %s \n",prelabel,value, postlabel); fi; fi; # End Function number 1 end; omniout_int := proc(iolevel,prelabel,prelen,value,vallen,postlabel) global glob_iolevel; if (glob_iolevel >= iolevel) then if vallen = 5 then printf("%-30s = %-32d %s\n",prelabel,value, postlabel); else printf("%-30s = %-32d %s \n",prelabel,value, postlabel); fi; fi; # End Function number 1 end; omniout_float_arr := proc(iolevel,prelabel,elemnt,prelen,value,vallen,postlabel) global glob_iolevel; if (glob_iolevel >= iolevel) then print(prelabel,"[",elemnt,"]",value, postlabel); fi; # End Function number 1 end; dump_series := proc(iolevel,dump_label,series_name, array_series,numb) global glob_iolevel; local i; if (glob_iolevel >= iolevel) then i := 1; while (i <= numb) do print(dump_label,series_name ,i,array_series[i]); i := i + 1; od; fi; # End Function number 1 end; dump_series_2 := proc(iolevel,dump_label,series_name2, array_series2,numb,subnum,array_x) global glob_iolevel; local i,sub,ts_term; if (glob_iolevel >= iolevel) then sub := 1; while (sub <= subnum) do i := 1; while (i <= numb) do print(dump_label,series_name2,sub,i,array_series2[sub,i]); od; sub := sub + 1; od; fi; # End Function number 1 end; cs_info := proc(iolevel,str) global glob_iolevel,glob_correct_start_flag,glob_h,glob_reached_optimal_h; if (glob_iolevel >= iolevel) then print("cs_info " , str , " glob_correct_start_flag = " , glob_correct_start_flag , "glob_h := " , glob_h , "glob_reached_optimal_h := " , glob_reached_optimal_h) fi; # End Function number 1 end; # Begin Function number 2 logitem_time := proc(fd,secs_in) global centuries_in_millinium, days_in_year, hours_in_day, min_in_hour, sec_in_min, years_in_century; local cent_int, centuries, days, days_int, hours, hours_int, millinium_int, milliniums, minutes, minutes_int, sec_in_millinium, sec_int, seconds, secs, years, years_int; secs := (secs_in); if (secs > 0.0) then # if number 1 sec_in_millinium := convfloat(sec_in_min * min_in_hour * hours_in_day * days_in_year * years_in_century * centuries_in_millinium); milliniums := convfloat(secs / sec_in_millinium); millinium_int := floor(milliniums); centuries := (milliniums - millinium_int)*centuries_in_millinium; cent_int := floor(centuries); years := (centuries - cent_int) * years_in_century; years_int := floor(years); days := (years - years_int) * days_in_year; days_int := floor(days); hours := (days - days_int) * hours_in_day; hours_int := floor(hours); minutes := (hours - hours_int) * min_in_hour; minutes_int := floor(minutes); seconds := (minutes - minutes_int) * sec_in_min; sec_int := floor(seconds); fprintf(fd,"