/* BEGIN OUTFILE1 */ // before write ccc cpp top matter #define true 1 #define false 0 #define ALWAYS 1 #define INFO 2 #define DEBUGL 3 #define DEBUGMASSIVE 4 #define glob_iolevel INFO #define MAX_TERMS 40 #include #include #include #include #include #define convfloat(x) ((double)x) #define arcsin(x) asin(x) #define arccos(x) acos(x) #define arctan(x) atan(x) #define float_abs(x) fabs(x) #define int_abs(x) abs(x) #define expt(x,y) pow(x,y) #define ln(x) log(x) #define Si(x) (0.0) #define Ci(x) (0.0) long elapsed_time_seconds(); int reached_interval(); int not_reached_end(double ,double ); double ats(int,double *,double *,int); double att(int,double *,double *,int); double factorial_1(int); double factorial_3(int,int); double comp_rad_from_ratio(double ,double ,int); double comp_ord_from_ratio(double ,double ,int); double comp_rad_from_three_terms(double ,double ,double ,int); double comp_ord_from_three_terms(double ,double ,double ,int); double comp_rad_from_six_terms(double ,double ,double ,double ,double ,double ,int); double comp_ord_from_six_terms(double ,double ,double ,double ,double ,double ,int); int int_trunc(double); void track_estimated_error(); void display_poles(); void display_pole_debug(int,int,double,double); double comp_expect_sec(double ,double ,double ,double ); double comp_percent(double ,double ,double ); void omniout_str(int iolevel,char *str); void omniout_str_noeol(int iolevel,char *str); void omniout_labstr(int iolevel,char *label,char *str); void omniout_float(int iolevel,char *prelabel,int prelen,double value,int vallen,char * postlabel); void omniout_int(int iolevel,char *prelabel,int prelen,int value,int vallen,char * postlabel); void dump_series(int iolevel,char *dump_label,char *series_name, double *array_series,int numb); void cs_info(int iolevel,char *str); void logitem_time(FILE *fd,double secs_in); double frac(double); void omniout_timestr(double secs_in); double array_const_m1[MAX_TERMS]; double estimated_needed_step_error(double x_start,double x_end,double glob_h,double est_answer); double test_suggested_h(); double est_size_answer(); double my_check_sign(double x0,double xf); double exact_soln_y(double ); double array_y_init[MAX_TERMS + 1]; double array_norms[MAX_TERMS + 1]; double array_fact_1[MAX_TERMS + 1]; double array_1st_rel_error[2 + 1]; double array_last_rel_error[2 + 1]; double array_est_rel_error[2 + 1]; double array_max_est_error[2 + 1]; int array_type_pole[2 + 1]; int array_type_real_pole[2 + 1]; int array_type_complex_pole[2 + 1]; int array_est_digits[2 + 1]; int array_good_digits[2 + 1]; double array_y[MAX_TERMS + 1]; double array_x[MAX_TERMS + 1]; double array_tmp0[MAX_TERMS + 1]; double array_tmp1_g[MAX_TERMS + 1]; double array_tmp1[MAX_TERMS + 1]; double array_tmp2[MAX_TERMS + 1]; double array_tmp3[MAX_TERMS + 1]; double array_m1[MAX_TERMS + 1]; double array_y_higher[2 + 1][MAX_TERMS + 1]; double array_y_higher_work[2 + 1][MAX_TERMS + 1]; double array_y_higher_work2[2 + 1][MAX_TERMS + 1]; double array_y_set_initial[2 + 1][MAX_TERMS + 1]; double array_given_rad_poles[2 + 1][3 + 1]; double array_given_ord_poles[2 + 1][3 + 1]; double array_rad_test_poles[2 + 1][4 + 1]; double array_ord_test_poles[2 + 1][4 + 1]; double array_fact_2[MAX_TERMS + 1][MAX_TERMS + 1]; // before generate constants // before generate globals definition /* Top Generate Globals Definition */ int MAX_UNCHANGED=10; double glob_prec=1.0e-16; int glob_est_digits=1; double glob_check_sign=1.0; double glob_desired_digits_correct=8.0; double glob_max_estimated_step_error=0.0; double glob_ratio_of_radius=0.1; double glob_percent_done=0.0; int glob_subiter_method=3; double glob_total_exp_sec=0.1; double glob_optimal_expect_sec=0.1; double glob_estimated_size_answer=100.0; unsigned char glob_html_log=true; int glob_good_digits=0; int glob_max_opt_iter=10; unsigned char glob_dump=false; unsigned char glob_djd_debug=true; unsigned char glob_display_flag=true; unsigned char glob_djd_debug2=true; int glob_h_reason=0; int glob_sec_in_minute=60; double glob_min_in_hour=60.0; double glob_hours_in_day=24.0; int glob_days_in_year=365; int glob_sec_in_hour=3600; int glob_sec_in_day=86400; int glob_sec_in_year=31536000; double glob_almost_1=0.9990; double glob_clock_sec=0.0; double glob_clock_start_sec=0.0; unsigned char glob_not_yet_finished=true; unsigned char glob_initial_pass=true; unsigned char glob_not_yet_start_msg=true; unsigned char glob_reached_optimal_h=false; unsigned char glob_optimal_done=false; double glob_disp_incr=0.1; double glob_h=0.1; double glob_diff_rc_fm=0.1; double glob_diff_rc_fmm1=0.1; double glob_diff_rc_fmm2=0.1; double glob_diff_ord_fm=0.1; double glob_diff_ord_fmm1=0.1; double glob_diff_ord_fmm2=0.1; double glob_six_term_ord_save=0.1; double glob_guess_error_rc=0.1; double glob_guess_error_ord=0.1; double glob_max_h=0.1; double glob_min_h=0.000001; int glob_type_given_pole=0; double glob_large_float=1.0e100; double glob_larger_float=1.1e100; double glob_least_given_sing=9.9e100; double glob_least_ratio_sing=9.9e100; double glob_least_3_sing=9.9e100; double glob_least_6_sing=9.9e100; double glob_last_good_h=0.1; unsigned char glob_look_poles=false; double glob_display_interval=0.0; double glob_next_display=0.0; unsigned char glob_dump_analytic=false; double glob_abserr=0.1e-10; double glob_relerr=0.1e-10; double glob_min_pole_est=0.1e+10; double glob_max_hours=0.0; int glob_max_iter=1000; double glob_max_rel_trunc_err=0.1e-10; double glob_max_trunc_err=0.1e-10; int glob_no_eqs=0; double glob_optimal_clock_start_sec=0.0; double glob_optimal_start=0.0; double glob_upper_ratio_limit=1.0001; double glob_lower_ratio_limit=0.9999; double glob_small_float=0.0; double glob_smallish_float=0.0; int glob_unchanged_h_cnt=0; unsigned char glob_warned=false; unsigned char glob_warned2=false; double glob_max_sec=10000.0; double glob_orig_start_sec=0.0; int glob_start=0; int glob_curr_iter_when_opt=0; int glob_current_iter=0; int glob_iter=0; double glob_normmax=0.0; double glob_max_minutes=0.0; /* Bottom Generate Globals Deninition */ // before generate const definition double array_const_1[MAX_TERMS + 1]; double array_const_0D0[MAX_TERMS + 1]; double array_const_2D0[MAX_TERMS + 1]; // before write_aux functions void display_poles() { double rad_given; if ((glob_type_given_pole == 1) || (glob_type_given_pole == 2)) { /* if number 1*/ rad_given = sqrt((array_x[1] - array_given_rad_poles[1][1]) * (array_x[1] - array_given_rad_poles[1][1]) + array_given_rad_poles[1][2] * array_given_rad_poles[1][2]); omniout_float(ALWAYS,"Radius of convergence (given) for eq 1 ",4,rad_given,4," "); omniout_float(ALWAYS,"Order of pole (given) ",4,array_given_ord_poles[1][1],4," "); if (rad_given < glob_least_given_sing) { /* if number 2*/ glob_least_given_sing = rad_given; }/* end if 2*/ ; } else if (glob_type_given_pole == 3) { /* if number 2*/ omniout_str(ALWAYS,"NO POLE (given) for Equation 1"); } else if (glob_type_given_pole == 5) { /* if number 3*/ omniout_str(ALWAYS,"SOME POLE (given) for Equation 1"); } else { omniout_str(ALWAYS,"NO INFO (given) for Equation 1"); }/* end if 3*/ ; if (array_rad_test_poles[1][1] < glob_large_float) { /* if number 3*/ omniout_float(ALWAYS,"Radius of convergence (ratio test) for eq 1 ",4,array_rad_test_poles[1][1],4," "); if (array_rad_test_poles[1][1]< glob_least_ratio_sing) { /* if number 4*/ glob_least_ratio_sing = array_rad_test_poles[1][1]; }/* end if 4*/ ; omniout_float(ALWAYS,"Order of pole (ratio test) ",4, array_ord_test_poles[1][1],4," "); } else { omniout_str(ALWAYS,"NO POLE (ratio test) for Equation 1"); }/* end if 3*/ ; if ((array_rad_test_poles[1][2] > 0.0) && (array_rad_test_poles[1][2] < glob_large_float)) { /* if number 3*/ omniout_float(ALWAYS,"Radius of convergence (three term test) for eq 1 ",4,array_rad_test_poles[1][2],4," "); if (array_rad_test_poles[1][2]< glob_least_3_sing) { /* if number 4*/ glob_least_3_sing = array_rad_test_poles[1][2]; }/* end if 4*/ ; omniout_float(ALWAYS,"Order of pole (three term test) ",4, array_ord_test_poles[1][2],4," "); } else { omniout_str(ALWAYS,"NO REAL POLE (three term test) for Equation 1"); }/* end if 3*/ ; if ((array_rad_test_poles[1][3] > 0.0) && (array_rad_test_poles[1][3] < glob_large_float)) { /* if number 3*/ omniout_float(ALWAYS,"Radius of convergence (six term test) for eq 1 ",4,array_rad_test_poles[1][3],4," "); if (array_rad_test_poles[1][3]< glob_least_6_sing) { /* if number 4*/ glob_least_6_sing = array_rad_test_poles[1][3]; }/* end if 4*/ ; omniout_float(ALWAYS,"Order of pole (six term test) ",4, array_ord_test_poles[1][3],4," "); } else { omniout_str(ALWAYS,"NO COMPLEX POLE (six term test) for Equation 1"); }/* end if 3*/ ; } /* End Function number1 */ double my_check_sign(double x0 ,double xf) { double ret; if (xf > x0) { /* if number 3*/ ret = 1.0; } else { ret = -1.0; }/* end if 3*/ ; return ret;; } /* End Function number1 */ double est_size_answer() { double min_size; min_size = glob_estimated_size_answer; if (float_abs(array_y[1]) < min_size) { /* if number 3*/ min_size = float_abs(array_y[1]); omniout_float(ALWAYS,"min_size",32,min_size,32,""); }/* end if 3*/ ; if (min_size < 1.0) { /* if number 3*/ min_size = 1.0; omniout_float(ALWAYS,"min_size",32,min_size,32,""); }/* end if 3*/ ; return min_size; } /* End Function number1 */ double test_suggested_h() { double max_estimated_step_error,hn_div_ho,hn_div_ho_2,hn_div_ho_3,est_tmp; int no_terms; max_estimated_step_error = 0.0; no_terms = MAX_TERMS; hn_div_ho = 0.5; hn_div_ho_2 = 0.25; hn_div_ho_3 = 0.125; omniout_float(ALWAYS,"hn_div_ho",32,hn_div_ho,32,""); omniout_float(ALWAYS,"hn_div_ho_2",32,hn_div_ho_2,32,""); omniout_float(ALWAYS,"hn_div_ho_3",32,hn_div_ho_3,32,""); est_tmp = float_abs(array_y[no_terms-3] + array_y[no_terms - 2] * hn_div_ho + array_y[no_terms - 1] * hn_div_ho_2 + array_y[no_terms] * hn_div_ho_3); if (est_tmp >= max_estimated_step_error) { /* if number 3*/ max_estimated_step_error = est_tmp; }/* end if 3*/ ; omniout_float(ALWAYS,"max_estimated_step_error",32,max_estimated_step_error,32,""); return max_estimated_step_error; } /* End Function number1 */ void track_estimated_error() { double hn_div_ho,hn_div_ho_2,hn_div_ho_3,est_tmp; int no_terms; no_terms = MAX_TERMS; hn_div_ho = 0.5; hn_div_ho_2 = 0.25; hn_div_ho_3 = 0.125; est_tmp = float_abs(array_y[no_terms-3]) + float_abs(array_y[no_terms - 2]) * hn_div_ho + float_abs(array_y[no_terms - 1]) * hn_div_ho_2 + float_abs(array_y[no_terms]) * hn_div_ho_3; if (glob_prec * float_abs(array_y[1]) > est_tmp) { /* if number 3*/ est_tmp = glob_prec * float_abs(array_y[1]); }/* end if 3*/ ; if (est_tmp >= array_max_est_error[1]) { /* if number 3*/ array_max_est_error[1] = est_tmp; }/* end if 3*/ ; } /* End Function number1 */ int reached_interval() { int ret; if ((glob_check_sign * array_x[1]) >= (glob_check_sign * glob_next_display)) { /* if number 3*/ ret = true; } else { ret = false; }/* end if 3*/ ; return(ret); } /* End Function number1 */ void display_alot(int iter) { double abserr, est_rel_err, analytic_val_y, ind_var, numeric_val, relerr; int term_no; /* TOP DISPLAY ALOT */ if (reached_interval()) { /* if number 3*/ if (iter >= 0) { /* if number 4*/ 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 = float_abs(numeric_val - analytic_val_y); omniout_float(ALWAYS,"y[1] (numeric) ",33,numeric_val,20," "); if (float_abs(analytic_val_y) > 0.000000000000000000000000000000001) { /* if number 5*/ relerr = abserr*100.0/float_abs(analytic_val_y); if (relerr > 0.000000000000000000000000000000001) { /* if number 6*/ glob_good_digits = -int_trunc(log10(relerr)) + 3; } else { glob_good_digits = 16; }/* end if 6*/ ; } else { relerr = -1.0 ; glob_good_digits = -16; }/* end if 5*/ ; if (float_abs(numeric_val) > 0.000000000000000000000000000000001) { /* if number 5*/ est_rel_err = array_max_est_error[1]*100.0 * sqrt(glob_iter)*19*MAX_TERMS/float_abs(numeric_val); if (est_rel_err > 0.000000000000000000000000000000001) { /* if number 6*/ glob_est_digits = -int_trunc(log10(est_rel_err)) + 2; } else { glob_est_digits = 16; }/* end if 6*/ ; } else { relerr = -1.0 ; glob_est_digits = -16; }/* end if 5*/ ; array_est_digits[1] = glob_est_digits; array_good_digits[1] = glob_good_digits; if (glob_iter == 1) { /* if number 5*/ array_1st_rel_error[1] = relerr; } else { array_last_rel_error[1] = relerr; }/* end if 5*/ ; array_est_rel_error[1] = est_rel_err; omniout_float(ALWAYS,"absolute error ",4,abserr,20," "); omniout_float(ALWAYS,"relative error ",4,relerr,20,"%"); omniout_int(INFO,"Desired digits ",32,glob_desired_digits_correct,4," "); omniout_int(INFO,"Estimated correct digits ",32,glob_est_digits,4," "); omniout_int(INFO,"Correct digits ",32,glob_good_digits,4," ") ; omniout_float(ALWAYS,"h ",4,glob_h,20," "); }/* end if 4*/ ; /* BOTTOM DISPLAY ALOT */ }/* end if 3*/ ; } /* End Function number1 */ void prog_report(double x_start,double x_end) { double clock_sec, opt_clock_sec, clock_sec1, expect_sec, left_sec, total_clock_sec; double percent_done; /* TOP PROGRESS REPORT */ clock_sec1 = elapsed_time_seconds(); total_clock_sec = (clock_sec1) - (glob_orig_start_sec); glob_clock_sec = (clock_sec1) - (glob_clock_start_sec); left_sec = (glob_max_sec) + (glob_orig_start_sec) - (clock_sec1); expect_sec = comp_expect_sec((x_end),(x_start),(array_x[1]) + (glob_h) ,( clock_sec1) - (glob_orig_start_sec)); opt_clock_sec = ( clock_sec1) - (glob_optimal_clock_start_sec); glob_optimal_expect_sec = comp_expect_sec((x_end),(x_start),(array_x[1]) +( glob_h) ,( opt_clock_sec)); glob_total_exp_sec = glob_optimal_expect_sec + total_clock_sec; percent_done = comp_percent((x_end),(x_start),(array_x[1]) + (glob_h)); glob_percent_done = percent_done; omniout_str_noeol(INFO,"Total Elapsed Time "); omniout_timestr((total_clock_sec)); omniout_str_noeol(INFO,"Elapsed Time(since restart) "); omniout_timestr((glob_clock_sec)); if ((percent_done) < (100.0)) { /* if number 3*/ omniout_str_noeol(INFO,"Expected Time Remaining "); omniout_timestr((expect_sec)); omniout_str_noeol(INFO,"Optimized Time Remaining "); omniout_timestr((glob_optimal_expect_sec)); omniout_str_noeol(INFO,"Expected Total Time "); omniout_timestr((glob_total_exp_sec)); }/* end if 3*/ ; omniout_str_noeol(INFO,"Time to Timeout "); omniout_timestr((left_sec)); omniout_float(INFO, "Percent Done ",33,percent_done,4,"%"); /* BOTTOM PROGRESS REPORT */ } /* End Function number1 */ void check_for_pole() { int cnt, m, n, found_sing, term, local_test, last_no; double dr1, dr2, ds1, ds2, hdrc, nr1, nr2, ord_no, term1, term2, term3, part1, part2, part3, part4, part5, part6, part7, part8, part9, part10, part11, part12, part13, part14, rad_c, rcs, rm0, rm1, rm2, rm3, rm4, h_new, ratio, tmp_rad,tmp_ord, tmp_ratio, prev_tmp_rad; /* TOP CHECK FOR POLE */ tmp_rad = glob_larger_float; prev_tmp_rad = glob_larger_float; tmp_ratio = glob_larger_float; rad_c = glob_larger_float; array_rad_test_poles[1][1] = glob_larger_float; array_ord_test_poles[1][1] = glob_larger_float; found_sing = 1; last_no = MAX_TERMS - 1 - 10; cnt = 0; while (last_no < MAX_TERMS-3 && found_sing == 1) { /* do number 1*/ tmp_rad = comp_rad_from_ratio(array_y_higher[1][last_no-1],array_y_higher[1][last_no],last_no); tmp_ratio = tmp_rad / prev_tmp_rad; if ((cnt > 0 ) && (tmp_ratio < glob_upper_ratio_limit) && (tmp_ratio > glob_lower_ratio_limit)) { /* if number 3*/ rad_c = tmp_rad; } else if (cnt == 0) { /* if number 4*/ rad_c = tmp_rad; } else if (cnt > 0) { /* if number 5*/ found_sing = 0; }/* end if 5*/ ; prev_tmp_rad = tmp_rad;; cnt = cnt + 1; last_no = last_no + 1; }/* end do number 1*/ ; if (found_sing == 1) { /* if number 5*/ if (rad_c < array_rad_test_poles[1][1]) { /* if number 6*/ array_rad_test_poles[1][1] = rad_c; last_no = last_no - 1; tmp_ord = comp_ord_from_ratio(array_y_higher[1][last_no-1],array_y_higher[1][last_no],last_no); array_rad_test_poles[1][1] = rad_c; array_ord_test_poles[1][1] = tmp_ord; }/* end if 6*/ ; }/* end if 5*/ ; /* BOTTOM general radius test1 */ tmp_rad = glob_larger_float; prev_tmp_rad = glob_larger_float; tmp_ratio = glob_larger_float; rad_c = glob_larger_float; array_rad_test_poles[1][2] = glob_larger_float; array_ord_test_poles[1][2] = glob_larger_float; found_sing = 1; last_no = MAX_TERMS - 1 - 10; cnt = 0; while (last_no < MAX_TERMS-4 && found_sing == 1) { /* do number 1*/ tmp_rad = comp_rad_from_three_terms(array_y_higher[1][last_no-2],array_y_higher[1][last_no-1],array_y_higher[1][last_no],last_no); tmp_ratio = tmp_rad / prev_tmp_rad; if ((cnt > 0 ) && (tmp_ratio < glob_upper_ratio_limit) && (tmp_ratio > glob_lower_ratio_limit)) { /* if number 5*/ rad_c = tmp_rad; } else if (cnt == 0) { /* if number 6*/ rad_c = tmp_rad; } else if (cnt > 0) { /* if number 7*/ found_sing = 0; }/* end if 7*/ ; prev_tmp_rad = tmp_rad;; cnt = cnt + 1; last_no = last_no + 1; }/* end do number 1*/ ; if (found_sing == 1) { /* if number 7*/ if (rad_c < array_rad_test_poles[1][2]) { /* if number 8*/ array_rad_test_poles[1][2] = rad_c; last_no = last_no - 1; tmp_ord = comp_ord_from_three_terms(array_y_higher[1][last_no-2],array_y_higher[1][last_no-1],array_y_higher[1][last_no],last_no); array_rad_test_poles[1][2] = rad_c; if (rad_c < glob_min_pole_est) { /* if number 9*/ glob_min_pole_est = rad_c; }/* end if 9*/ ; array_ord_test_poles[1][2] = tmp_ord; }/* end if 8*/ ; }/* end if 7*/ ; /* BOTTOM general radius test1 */ tmp_rad = glob_larger_float; prev_tmp_rad = glob_larger_float; tmp_ratio = glob_larger_float; rad_c = glob_larger_float; array_rad_test_poles[1][3] = glob_larger_float; array_ord_test_poles[1][3] = glob_larger_float; found_sing = 1; last_no = MAX_TERMS - 1 - 10; cnt = 0; while (last_no < MAX_TERMS-7 && found_sing == 1) { /* do number 1*/ tmp_rad = comp_rad_from_six_terms(array_y_higher[1][last_no-5],array_y_higher[1][last_no-4],array_y_higher[1][last_no-3],array_y_higher[1][last_no-2],array_y_higher[1][last_no-1],array_y_higher[1][last_no],last_no); tmp_ratio = tmp_rad / prev_tmp_rad; if ((cnt > 0 ) && (tmp_ratio < glob_upper_ratio_limit) && (tmp_ratio > glob_lower_ratio_limit)) { /* if number 7*/ rad_c = tmp_rad; } else if (cnt == 0) { /* if number 8*/ rad_c = tmp_rad; } else if (cnt > 0) { /* if number 9*/ found_sing = 0; }/* end if 9*/ ; prev_tmp_rad = tmp_rad;; cnt = cnt + 1; last_no = last_no + 1; }/* end do number 1*/ ; if (found_sing == 1) { /* if number 9*/ if (rad_c < array_rad_test_poles[1][3]) { /* if number 10*/ array_rad_test_poles[1][3] = rad_c; last_no = last_no - 1; tmp_ord = comp_ord_from_six_terms(array_y_higher[1][last_no-5],array_y_higher[1][last_no-4],array_y_higher[1][last_no-3],array_y_higher[1][last_no-2],array_y_higher[1][last_no-1],array_y_higher[1][last_no],last_no); array_rad_test_poles[1][3] = rad_c; if (rad_c < glob_min_pole_est) { /* if number 11*/ glob_min_pole_est = rad_c; }/* end if 11*/ ; array_ord_test_poles[1][3] = tmp_ord; }/* end if 10*/ ; }/* end if 9*/ ; /* BOTTOM general radius test1 */ /* START ADJUST ALL SERIES */ if (float_abs(glob_min_pole_est) * glob_ratio_of_radius < float_abs(glob_h)) { /* if number 9*/ h_new = glob_check_sign * glob_min_pole_est * glob_ratio_of_radius; omniout_str(ALWAYS,"SETTING H FOR POLE"); glob_h_reason = 6; if (glob_check_sign * glob_min_h > glob_check_sign * h_new) { /* if number 10*/ omniout_str(ALWAYS,"SETTING H FOR MIN H"); h_new = glob_min_h; glob_h_reason = 5; }/* end if 10*/ ; term = 1; ratio = 1.0; while (term <= MAX_TERMS) { /* do number 1*/ array_y[term] = array_y[term]* ratio; array_y_higher[1][term] = array_y_higher[1][term]* ratio; array_x[term] = array_x[term]* ratio; ratio = ratio * h_new / float_abs(glob_h); term = term + 1; }/* end do number 1*/ ; glob_h = h_new; }/* end if 9*/ ; /* BOTTOM ADJUST ALL SERIES */ ; if (reached_interval()) { /* if number 9*/ display_poles(); }/* end if 9*/ } /* End Function number1 */ void atomall() { int kkk, order_d, term, adj2, adj3; double temporary,temp,temp2; /* TOP ATOMALL */ // before write ccc cpp library user def /* END OUTFILE1 */ /* BEGIN OUTFILE2 */ /* END OUTFILE2 */ /* BEGIN ATOMHDR1 */ /* emit pre sin 1 $eq_no = 1 */ array_tmp1[1] = sin(array_x[1]); array_tmp1_g[1] = cos(array_x[1]); /* emit pre mult CONST FULL $eq_no = 1 i = 1 */ array_tmp2[1] = array_const_2D0[1] * array_tmp1[1]; /* emit pre add CONST FULL $eq_no = 1 i = 1 */ array_tmp3[1] = array_const_0D0[1] + array_tmp2[1]; /* emit pre assign xxx $eq_no = 1 i = 1 $min_hdrs = 5 */ if ( ! array_y_set_initial[1][2]) { /* if number 1*/ if (1 <= MAX_TERMS) { /* if number 2*/ temporary = array_tmp3[1] * expt(glob_h , (1)) * factorial_3(0,1); if (2 <= MAX_TERMS) { /* if number 3*/ array_y[2] = temporary; array_y_higher[1][2] = temporary; }/* end if 3*/ ; temporary = temporary / glob_h * (1.0); array_y_higher[2][1] = temporary; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 2; /* END ATOMHDR1 */ /* BEGIN ATOMHDR2 */ /* emit pre sin ID_LINEAR iii = 2 $eq_no = 1 */ array_tmp1[2] = array_tmp1_g[1] * array_x[2] / 1; array_tmp1_g[2] = -array_tmp1[1] * array_x[2] / 1; /* emit pre mult CONST FULL $eq_no = 1 i = 2 */ array_tmp2[2] = array_const_2D0[1] * array_tmp1[2]; /* emit pre add CONST FULL $eq_no = 1 i = 2 */ array_tmp3[2] = array_tmp2[2]; /* emit pre assign xxx $eq_no = 1 i = 2 $min_hdrs = 5 */ if ( ! array_y_set_initial[1][3]) { /* if number 1*/ if (2 <= MAX_TERMS) { /* if number 2*/ temporary = array_tmp3[2] * expt(glob_h , (1)) * factorial_3(1,2); if (3 <= MAX_TERMS) { /* if number 3*/ array_y[3] = temporary; array_y_higher[1][3] = temporary; }/* end if 3*/ ; temporary = temporary / glob_h * (2.0); array_y_higher[2][2] = temporary; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 3; /* END ATOMHDR2 */ /* BEGIN ATOMHDR3 */ /* emit pre sin ID_LINEAR iii = 3 $eq_no = 1 */ array_tmp1[3] = array_tmp1_g[2] * array_x[2] / 2; array_tmp1_g[3] = -array_tmp1[2] * array_x[2] / 2; /* emit pre mult CONST FULL $eq_no = 1 i = 3 */ array_tmp2[3] = array_const_2D0[1] * array_tmp1[3]; /* emit pre add CONST FULL $eq_no = 1 i = 3 */ array_tmp3[3] = array_tmp2[3]; /* emit pre assign xxx $eq_no = 1 i = 3 $min_hdrs = 5 */ if ( ! array_y_set_initial[1][4]) { /* if number 1*/ if (3 <= MAX_TERMS) { /* if number 2*/ temporary = array_tmp3[3] * expt(glob_h , (1)) * factorial_3(2,3); if (4 <= MAX_TERMS) { /* if number 3*/ array_y[4] = temporary; array_y_higher[1][4] = temporary; }/* end if 3*/ ; temporary = temporary / glob_h * (3.0); array_y_higher[2][3] = temporary; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 4; /* END ATOMHDR3 */ /* BEGIN ATOMHDR4 */ /* emit pre sin ID_LINEAR iii = 4 $eq_no = 1 */ array_tmp1[4] = array_tmp1_g[3] * array_x[2] / 3; array_tmp1_g[4] = -array_tmp1[3] * array_x[2] / 3; /* emit pre mult CONST FULL $eq_no = 1 i = 4 */ array_tmp2[4] = array_const_2D0[1] * array_tmp1[4]; /* emit pre add CONST FULL $eq_no = 1 i = 4 */ array_tmp3[4] = array_tmp2[4]; /* emit pre assign xxx $eq_no = 1 i = 4 $min_hdrs = 5 */ if ( ! array_y_set_initial[1][5]) { /* if number 1*/ if (4 <= MAX_TERMS) { /* if number 2*/ temporary = array_tmp3[4] * expt(glob_h , (1)) * factorial_3(3,4); if (5 <= MAX_TERMS) { /* if number 3*/ array_y[5] = temporary; array_y_higher[1][5] = temporary; }/* end if 3*/ ; temporary = temporary / glob_h * (4.0); array_y_higher[2][4] = temporary; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 5; /* END ATOMHDR4 */ /* BEGIN ATOMHDR5 */ /* emit pre sin ID_LINEAR iii = 5 $eq_no = 1 */ array_tmp1[5] = array_tmp1_g[4] * array_x[2] / 4; array_tmp1_g[5] = -array_tmp1[4] * array_x[2] / 4; /* emit pre mult CONST FULL $eq_no = 1 i = 5 */ array_tmp2[5] = array_const_2D0[1] * array_tmp1[5]; /* emit pre add CONST FULL $eq_no = 1 i = 5 */ array_tmp3[5] = array_tmp2[5]; /* emit pre assign xxx $eq_no = 1 i = 5 $min_hdrs = 5 */ if ( ! array_y_set_initial[1][6]) { /* if number 1*/ if (5 <= MAX_TERMS) { /* if number 2*/ temporary = array_tmp3[5] * expt(glob_h , (1)) * factorial_3(4,5); if (6 <= MAX_TERMS) { /* if number 3*/ array_y[6] = temporary; array_y_higher[1][6] = temporary; }/* end if 3*/ ; temporary = temporary / glob_h * (5.0); array_y_higher[2][5] = temporary; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 6; /* END ATOMHDR5 */ /* BEGIN OUTFILE3 */ /* Top Atomall While Loop-- outfile3 */ while (kkk <= MAX_TERMS) { /* do number 1*/ /* END OUTFILE3 */ /* BEGIN OUTFILE4 */ /* emit sin LINEAR $eq_no = 1 */ array_tmp1[kkk] = array_tmp1_g[kkk - 1] * array_x[2] / (kkk - 1); array_tmp1_g[kkk] = -array_tmp1[kkk - 1] * array_x[2] / (kkk - 1); /* emit mult CONST FULL $eq_no = 1 i = 1 */ array_tmp2[kkk] = array_const_2D0[1] * array_tmp1[kkk]; /* emit NOT FULL - FULL add $eq_no = 1 */ array_tmp3[kkk] = array_tmp2[kkk]; /* emit assign $eq_no = 1 */ order_d = 1; if (kkk + order_d <= MAX_TERMS) { /* if number 1*/ if ( ! array_y_set_initial[1][kkk + order_d]) { /* if number 2*/ temporary = array_tmp3[kkk] * expt(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 = kkk + order_d - 1; adj3 = 2; while ((term >= 1) && (term <= MAX_TERMS) && (adj3 < order_d + 1)) { /* do number 1*/ if (adj3 <= order_d + 1) { /* if number 3*/ if (adj2 > 0) { /* if number 4*/ temporary = temporary / glob_h * (adj2); } else { temporary = temporary; }/* end if 4*/ ; array_y_higher[adj3][term] = temporary; }/* end if 3*/ ; term = term - 1; adj2 = adj2 - 1; adj3 = adj3 + 1; }/* end do number 1*/ }/* end if 2*/ }/* end if 1*/ ; kkk = kkk + 1; }/* end do number 1*/ ; /* BOTTOM ATOMALL */ /* END OUTFILE4 */ /* BEGIN OUTFILE5 */ /* BOTTOM ATOMALL ??? */ } /* End Function number1 */ /* END OUTFILE5 */ // before generate arrays /* BEGIN ATS LIBRARY BLOCK */ void omniout_str(int iolevel,char *str) { if (glob_iolevel >= iolevel) { printf("%s\n",str); } } /* End Function number1 */ void omniout_str_noeol(int iolevel,char *str) { if (glob_iolevel >= iolevel) { printf("%s",str); } } /* End Function number1 */ void omniout_labstr(int iolevel,char *label,char *str) { if (glob_iolevel >= iolevel) { printf("%s = %s\n",label,str); } } /* End Function number1 */ void omniout_float(int iolevel,char *prelabel,int prelen,double value,int vallen,char *postlabel) { if (glob_iolevel >= iolevel) { if (vallen == 4) { printf("%-30s = %-42.4g %s \n",prelabel,value, postlabel); } else { printf("%-30s = %-42.16g %s \n",prelabel,value, postlabel); } } } /* End Function number1 */ void omniout_int(int iolevel,char *prelabel,int prelen,int value,int vallen,char *postlabel) { if (glob_iolevel >= iolevel) { if (vallen == 5) { printf("%-30s = %-32d %s\n",prelabel,value, postlabel); } else { printf("%-30s = %-32d %s \n",prelabel,value, postlabel); } } } /* End Function number1 */ void omniout_float_arr(int iolevel,char *prelabel,int elemnt,int prelen,double *value,int vallen,char *postlabel) { if (glob_iolevel >= iolevel) { printf("%s = [ %d ] %g %s \n", prelabel,elemnt,value, postlabel); } } /* End Function number1 */ void logitem_time(FILE *fd,double secs_in) { int days_int, hours_int, minutes_int, sec_int, years_int,sec_temp; double sec_dbl; fprintf(fd,""); if (secs_in >= 0) { /* if number 1*/ years_int = (int_trunc(secs_in) / glob_sec_in_year); sec_temp = (int_trunc(secs_in) % glob_sec_in_year); days_int = (sec_temp / glob_sec_in_day) ; sec_temp = (sec_temp % glob_sec_in_day) ; hours_int = (sec_temp / glob_sec_in_hour); sec_temp = (sec_temp % glob_sec_in_hour); minutes_int = (sec_temp / glob_sec_in_minute); sec_int = (sec_temp % glob_sec_in_minute); sec_dbl = secs_in - (double)(years_int * glob_sec_in_year + days_int * glob_sec_in_day + hours_int * glob_sec_in_hour + minutes_int * glob_sec_in_minute); if (years_int > 0) { /* if number 2*/ fprintf(fd,"%d Years %d Days %d Hours %d Minutes %3.1f Seconds",years_int,days_int,hours_int,minutes_int,sec_dbl); } else if (days_int > 0) { /* if number 3*/ fprintf(fd,"%d Days %d Hours %d Minutes %3.1f Seconds",days_int,hours_int,minutes_int,sec_dbl); } else if (hours_int > 0) { /* if number 4*/ fprintf(fd,"%d Hours %d Minutes %3.1f Seconds",hours_int,minutes_int,sec_dbl); } else if (minutes_int > 0) { /* if number 5*/ fprintf(fd,"%d Minutes %3.1f Seconds",minutes_int,sec_dbl); } else { fprintf(fd,"%3.1f Seconds",sec_dbl); }/* end if 5*/ } else { fprintf(fd,"0.0 Seconds"); }/* end if 4*/ fprintf(fd,""); } /* End Function number1 */ void omniout_timestr(double secs_in) { int days_int, hours_int, minutes_int, sec_int, years_int,sec_temp; double sec_dbl; if (secs_in >= 0) { /* if number 4*/ years_int = (int_trunc(secs_in) / glob_sec_in_year); sec_temp = (int_trunc(secs_in) % glob_sec_in_year); days_int = (sec_temp / glob_sec_in_day) ; sec_temp = (sec_temp % glob_sec_in_day) ; hours_int = (sec_temp / glob_sec_in_hour); sec_temp = (sec_temp % glob_sec_in_hour); minutes_int = (sec_temp / glob_sec_in_minute); sec_int = (sec_temp % glob_sec_in_minute); sec_dbl = secs_in - (double)(years_int * glob_sec_in_year + days_int * glob_sec_in_day + hours_int * glob_sec_in_hour + minutes_int * glob_sec_in_minute); if (years_int > 0) { /* if number 5*/ printf(" = %d Years %d Days %d Hours %d Minutes %3.1f Seconds\n",years_int,days_int,hours_int,minutes_int,sec_dbl); } else if (days_int > 0) { /* if number 6*/ printf(" = %d Days %d Hours %d Minutes %3.1f Seconds\n",days_int,hours_int,minutes_int,sec_dbl); } else if (hours_int > 0) { /* if number 7*/ printf(" = %d Hours %d Minutes %3.1f Seconds\n",hours_int,minutes_int,sec_dbl); } else if (minutes_int > 0) { /* if number 8*/ printf(" = %d Minutes %3.1f Seconds\n",minutes_int,sec_dbl); } else { printf(" = %3.1f Seconds\n",sec_dbl); }/* end if 8*/ } else { printf(" 0.0 Seconds\n"); }/* end if 7*/ } /* End Function number1 */ double zero_ats_ar(double *arr_a) { int iii; iii = 1; while (iii <= MAX_TERMS) { /* do number 1*/ arr_a [iii] = 0.0; iii = iii + 1; }/* end do number 1*/ } /* End Function number1 */ double ats(int mmm_ats,double *arr_a,double *arr_b,int jjj_ats) { int iii_ats, lll_ats, ma_ats; double ret_ats; ret_ats = 0.0; if (jjj_ats <= mmm_ats) { /* if number 7*/ ma_ats = mmm_ats + 1; iii_ats = jjj_ats; while (iii_ats <= mmm_ats) { /* do number 1*/ lll_ats = ma_ats - iii_ats; if ((lll_ats <= MAX_TERMS && (iii_ats <= MAX_TERMS) )) { /* if number 8*/ ret_ats = ret_ats + arr_a[iii_ats]*arr_b[lll_ats]; }/* end if 8*/ ; iii_ats = iii_ats + 1; }/* end do number 1*/ }/* end if 7*/ ; return ret_ats; } /* End Function number1 */ double att(int mmm_att,double *arr_aa,double *arr_bb,int jjj_att) { int al_att, iii_att,lll_att, ma_att; double ret_att; ret_att = 0.0; if (jjj_att < mmm_att) { /* if number 7*/ ma_att = mmm_att + 2; iii_att = jjj_att; while ((iii_att < mmm_att) && (iii_att <= MAX_TERMS) ) { /* do number 1*/ lll_att = ma_att - iii_att; al_att = (lll_att - 1); if ((lll_att <= MAX_TERMS && (iii_att <= MAX_TERMS) )) { /* if number 8*/ ret_att = ret_att + arr_aa[iii_att]*arr_bb[lll_att]* (al_att); }/* end if 8*/ ; iii_att = iii_att + 1; }/* end do number 1*/ ; ret_att = ret_att / (mmm_att) ; }/* end if 7*/ ; return ret_att; } /* End Function number1 */ void logditto(FILE * file) { fprintf(file,""); fprintf(file,"ditto"); fprintf(file,""); } /* End Function number1 */ void logitem_integer(FILE * file,int n) { fprintf(file,""); fprintf(file,"%d",n); fprintf(file,""); } /* End Function number1 */ void logitem_str(FILE * file,char * str) { fprintf(file,""); fprintf(file,str); fprintf(file,""); } /* End Function number1 */ void logitem_good_digits(FILE * file,double rel_error) { int good_digits; fprintf(file,""); if (rel_error != -1.0) { /* if number 7*/ if (rel_error > 0.000000000000000000000000000000001) { /* if number 8*/ good_digits = 3-int_trunc(log10(rel_error)); fprintf(file,"%d",good_digits); } else { good_digits = 16; fprintf(file,"%d",good_digits); }/* end if 8*/ ; } else { fprintf(file,"Unknown"); }/* end if 7*/ ; fprintf(file,""); } /* End Function number1 */ void log_revs(FILE * file,char * revs) { fprintf(file,revs); } /* End Function number1 */ void logitem_float(FILE * file,double x) { fprintf(file,""); fprintf(file,"%g",x); fprintf(file,""); } /* End Function number1 */ void logitem_h_reason(FILE * file) { fprintf(file,""); if (glob_h_reason == 1) { /* if number 7*/ fprintf(file,"Max H"); } else if (glob_h_reason == 2) { /* if number 8*/ fprintf(file,"Display Interval"); } else if (glob_h_reason == 3) { /* if number 9*/ fprintf(file,"Optimal"); } else if (glob_h_reason == 4) { /* if number 10*/ fprintf(file,"Pole Accuracy"); } else if (glob_h_reason == 5) { /* if number 11*/ fprintf(file,"Min H (Pole)"); } else if (glob_h_reason == 6) { /* if number 12*/ fprintf(file,"Pole"); } else { fprintf(file,"Impossible"); }/* end if 12*/ fprintf(file,""); } /* End Function number1 */ void logstart(FILE * file) { fprintf(file,""); } /* End Function number1 */ void logend(FILE * file) { fprintf(file,"\n"); } /* End Function number1 */ void chk_data() { int errflag; errflag = false; if (glob_max_iter < 2) { /* if number 12*/ omniout_str(ALWAYS,"Illegal max_iter"); errflag = true; }/* end if 12*/ ; if (errflag) { /* if number 12*/ }/* end if 12*/ } /* End Function number1 */ double comp_expect_sec(double t_end2,double t_start2,double t2,double clock_sec2) { double ms2, rrr, sec_left, sub1, sub2; ; ms2 = clock_sec2; sub1 = (t_end2-t_start2); sub2 = (t2-t_start2); if (sub1 == 0.0) { /* if number 12*/ sec_left = 0.0; } else { if (sub2 > 0.0) { /* if number 13*/ rrr = (sub1/sub2); sec_left = rrr * ms2 - ms2; } else { sec_left = 0.0; }/* end if 13*/ }/* end if 12*/ ; return sec_left; } /* End Function number1 */ double comp_percent(double t_end2,double t_start2,double t2) { double rrr, sub1, sub2; sub1 = (t_end2-t_start2); sub2 = (t2-t_start2); if (sub2 > glob_small_float) { /* if number 12*/ rrr = (100.0*sub2)/sub1; } else { rrr = 0.0; }/* end if 12*/ ; return rrr; } /* End Function number1 */ double comp_rad_from_ratio(double term1,double term2,int last_no) { /* TOP TWO TERM RADIUS ANALYSIS */ double ret; if (term2 > 0.0) { /* if number 12*/ ret = float_abs(term1 * glob_h / term2); } else { ret = glob_larger_float; }/* end if 12*/ ; return ret; /* BOTTOM TWO TERM RADIUS ANALYSIS */ } /* End Function number1 */ double comp_ord_from_ratio(double term1,double term2,int last_no) { /* TOP TWO TERM ORDER ANALYSIS */ double ret; if (term2 > 0.0) { /* if number 12*/ ret = 1.0 + float_abs(term2) * convfloat(last_no) * ln(float_abs(term1 * glob_h / term2))/ln(convfloat(last_no)); } else { ret = glob_larger_float; }/* end if 12*/ ; return ret; /* BOTTOM TWO TERM ORDER ANALYSIS */ } /* End Function number1 */ double comp_rad_from_three_terms(double term1,double term2,double term3,int last_no) { /* TOP THREE TERM RADIUS ANALYSIS */ double ret,temp; temp = float_abs(term2*term2*convfloat(last_no)-2.0*term2*term2-term1*term3*convfloat(last_no)+term1*term3); if (temp > 0.0) { /* if number 12*/ ret = float_abs((term2*glob_h*term1)/(temp)); } else { ret = glob_larger_float; }/* end if 12*/ ; return ret; /* BOTTOM THREE TERM RADIUS ANALYSIS */ } /* End Function number1 */ double comp_ord_from_three_terms(double term1,double term2,double term3,int last_no) { /* TOP THREE TERM ORDER ANALYSIS */ double ret; ret = float_abs((4.0*term1*term3*convfloat(last_no)-3.0*term1*term3-4.0*term2*term2*convfloat(last_no)+4.0*term2*term2+term2*term2*convfloat(last_no*last_no)-term1*term3*convfloat(last_no*last_no))/(term2*term2*convfloat(last_no)-2.0*term2*term2-term1*term3*convfloat(last_no)+term1*term3)); return ret; /* TOP THREE TERM ORDER ANALYSIS */ } /* End Function number1 */ double comp_rad_from_six_terms(double term1,double term2,double term3,double term4,double term5,double term6,int last_no) { /* TOP SIX TERM RADIUS ANALYSIS */ double ret,rm0,rm1,rm2,rm3,rm4,nr1,nr2,dr1,dr2,ds2,rad_c,ord_no,ds1,rcs; if ((term5 != 0.0) && (term4 != 0.0) && (term3 != 0.0) && (term2 != 0.0) && (term1 != 0.0)) { /* if number 12*/ rm0 = term6/term5; rm1 = term5/term4; rm2 = term4/term3; rm3 = term3/term2; rm4 = term2/term1; nr1 = convfloat(last_no-1)*rm0 - 2.0*convfloat(last_no-2)*rm1 + convfloat(last_no-3)*rm2; nr2 = convfloat(last_no-2)*rm1 - 2.0*convfloat(last_no-3)*rm2 + convfloat(last_no-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 ((float_abs(nr1 * dr2 - nr2 * dr1) == 0.0) || (float_abs(dr1) == 0.0)) { /* if number 13*/ rad_c = glob_larger_float; ord_no = glob_larger_float; } else { if (float_abs(nr1*dr2 - nr2 * dr1) != 0.0) { /* if number 14*/ 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(last_no)/2.0; if (float_abs(rcs) != 0.0) { /* if number 15*/ if (rcs > 0.0) { /* if number 16*/ rad_c = sqrt(rcs) * float_abs(glob_h); } else { rad_c = glob_larger_float; ord_no = glob_larger_float; }/* end if 16*/ } else { rad_c = glob_larger_float; ord_no = glob_larger_float; }/* end if 15*/ } else { rad_c = glob_larger_float; ord_no = glob_larger_float; }/* end if 14*/ }/* end if 13*/ } else { rad_c = glob_larger_float; ord_no = glob_larger_float; }/* end if 12*/ ; glob_six_term_ord_save = ord_no; return rad_c; /* BOTTOM SIX TERM RADIUS ANALYSIS */ } /* End Function number1 */ double comp_ord_from_six_terms(double term1,double term2,double term3,double term4,double term5,double term6,int last_no) { double ret; /* TOP SIX TERM ORDER ANALYSIS */ /* TOP SAVED FROM SIX TERM RADIUS ANALYSIS */ return glob_six_term_ord_save; /* BOTTOM SIX TERM ORDER ANALYSIS */ } /* End Function number1 */ double factorial_2(int nnn) { double ret; if (nnn > 0) { /* if number 12*/ return nnn * factorial_2(nnn - 1); } else { return 1.0; }/* end if 12*/ } /* End Function number1 */ double factorial_1(int nnn) { double ret; if (nnn <= MAX_TERMS) { /* if number 12*/ if (array_fact_1[nnn] == 0) { /* if number 13*/ ret = factorial_2(nnn); array_fact_1[nnn] = ret; } else { ret = array_fact_1[nnn]; }/* end if 13*/ ; } else { ret = factorial_2(nnn); }/* end if 12*/ ; return ret; } /* End Function number1 */ double factorial_3(int mmm,int nnn) { double ret; if ((nnn <= MAX_TERMS) && (mmm <= MAX_TERMS)) { /* if number 12*/ if (array_fact_2[mmm][nnn] == 0) { /* if number 13*/ ret = factorial_1(mmm)/factorial_1(nnn); array_fact_2[mmm][nnn] = ret; } else { ret = array_fact_2[mmm][nnn]; }/* end if 13*/ ; } else { ret = factorial_2(mmm)/factorial_2(nnn); }/* end if 12*/ ; return ret; } /* End Function number1 */ double frac(double in) { double out; out = in - trunc(in); return(out); } /* End Function number1 */ int int_trunc(double in) { int out; out = (int)trunc(in); return(out); } /* End Function number1 */ long elapsed_time_seconds() { struct timeval t; struct timezone tz; gettimeofday(&t,&tz); return (t.tv_sec); } /* End Function number1 */ double estimated_needed_step_error(double x_start,double x_end,double estimated_h,double estimated_answer) { double desired_abs_gbl_error,range,estimated_steps,step_error; omniout_float(ALWAYS,"glob_desired_digits_correct",32,glob_desired_digits_correct,32,""); desired_abs_gbl_error = expt(10.0, -glob_desired_digits_correct) * float_abs(estimated_answer); omniout_float(ALWAYS,"estimated_h",32,estimated_h,32,""); omniout_float(ALWAYS,"estimated_answer",32,estimated_answer,32,""); omniout_float(ALWAYS,"desired_abs_gbl_error",32,desired_abs_gbl_error,32,""); range = (x_end - x_start); omniout_float(ALWAYS,"range",32,range,32,""); estimated_steps = range / estimated_h; omniout_float(ALWAYS,"estimated_steps",32,estimated_steps,32,""); step_error = (float_abs(desired_abs_gbl_error /sqrt( estimated_steps)/MAX_TERMS)); omniout_float(ALWAYS,"step_error",32,step_error,32,""); return (step_error);; } /* End Function number1 */ /* END ATS LIBRARY BLOCK */ /* BEGIN USER FUNCTION BLOCK */ /* BEGIN BLOCK 3 */ /* BEGIN USER DEF BLOCK */ double exact_soln_y (double x) { return(-cos(x) * 2.0); } /* END USER DEF BLOCK */ /* END BLOCK 3 */ /* END USER FUNCTION BLOCK */ int main() { /* BEGIN OUTFIEMAIN */ double d1,d2,d3,d4,est_err_2,Digits; int niii,done_once,it,opt_iter, max_terms; int term,ord,order_diff,term_no,iiif,jjjf,subiter; FILE *html_log_file; int rows,r_order,sub_iter,calc_term,iii,current_iter,repeat_it,found_h; double temp_sum, est_needed_step_err,estimated_step_error,min_value,est_answer,last_min_pole_est; double x_start;double x_end; double tmp; // before first input block /* BEGIN FIRST INPUT BLOCK */ /* BEGIN BLOCK 1 */ /* BEGIN FIRST INPUT BLOCK */ Digits=32; max_terms=40; /* END BLOCK 1 */ /* END FIRST INPUT BLOCK */ /* START OF INITS AFTER INPUT BLOCK */ glob_html_log = true; /* END OF INITS AFTER INPUT BLOCK */ // before arrays initialized term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_y_init[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_norms[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_fact_1[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= 2) { /* do number 1*/ array_1st_rel_error[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= 2) { /* do number 1*/ array_last_rel_error[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= 2) { /* do number 1*/ array_est_rel_error[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= 2) { /* do number 1*/ array_max_est_error[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= 2) { /* do number 1*/ array_type_pole[term] = 0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= 2) { /* do number 1*/ array_type_real_pole[term] = 0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= 2) { /* do number 1*/ array_type_complex_pole[term] = 0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= 2) { /* do number 1*/ array_est_digits[term] = 0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= 2) { /* do number 1*/ array_good_digits[term] = 0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_y[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_x[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_tmp0[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_tmp1_g[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_tmp1[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_tmp2[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_tmp3[term] = 0.0; term = term + 1; }/* end do number 1*/ ; term = 1; while (term <= MAX_TERMS) { /* do number 1*/ array_m1[term] = 0.0; term = term + 1; }/* end do number 1*/ ; ord = 1; while (ord <=2) { /* do number 1*/ term = 1; while (term <= MAX_TERMS) { /* do number 2*/ array_y_higher[ord][term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = ord + 1; }/* end do number 1*/ ; ord = 1; while (ord <=2) { /* do number 1*/ term = 1; while (term <= MAX_TERMS) { /* do number 2*/ array_y_higher_work[ord][term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = ord + 1; }/* end do number 1*/ ; ord = 1; while (ord <=2) { /* do number 1*/ term = 1; while (term <= MAX_TERMS) { /* do number 2*/ array_y_higher_work2[ord][term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = ord + 1; }/* end do number 1*/ ; ord = 1; while (ord <=2) { /* do number 1*/ term = 1; while (term <= MAX_TERMS) { /* do number 2*/ array_y_set_initial[ord][term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = ord + 1; }/* end do number 1*/ ; ord = 1; while (ord <=2) { /* do number 1*/ term = 1; while (term <= 3) { /* do number 2*/ array_given_rad_poles[ord][term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = ord + 1; }/* end do number 1*/ ; ord = 1; while (ord <=2) { /* do number 1*/ term = 1; while (term <= 3) { /* do number 2*/ array_given_ord_poles[ord][term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = ord + 1; }/* end do number 1*/ ; ord = 1; while (ord <=2) { /* do number 1*/ term = 1; while (term <= 4) { /* do number 2*/ array_rad_test_poles[ord][term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = ord + 1; }/* end do number 1*/ ; ord = 1; while (ord <=2) { /* do number 1*/ term = 1; while (term <= 4) { /* do number 2*/ array_ord_test_poles[ord][term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = ord + 1; }/* end do number 1*/ ; ord = 1; while (ord <=MAX_TERMS) { /* do number 1*/ term = 1; while (term <= MAX_TERMS) { /* do number 2*/ array_fact_2[ord][term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = ord + 1; }/* end do number 1*/ ; // before symbols initialized /* BEGIN SYMBOLS INITIALIZATED */ zero_ats_ar(array_y); zero_ats_ar(array_x); zero_ats_ar(array_tmp0); zero_ats_ar(array_tmp1_g); zero_ats_ar(array_tmp1); zero_ats_ar(array_tmp2); zero_ats_ar(array_tmp3); zero_ats_ar(array_m1); zero_ats_ar(array_const_1); array_const_1[1] = 1; zero_ats_ar(array_const_0D0); array_const_0D0[1] = 0.0; zero_ats_ar(array_const_2D0); array_const_2D0[1] = 2.0; zero_ats_ar(array_m1); array_m1[1] = -1.0; /* END SYMBOLS INITIALIZATED */ MAX_UNCHANGED = 10; glob_prec = 1.0e-16; glob_est_digits = 1; glob_check_sign = 1.0; glob_desired_digits_correct = 8.0; glob_max_estimated_step_error = 0.0; glob_ratio_of_radius = 0.1; glob_percent_done = 0.0; glob_subiter_method = 3; glob_total_exp_sec = 0.1; glob_optimal_expect_sec = 0.1; glob_estimated_size_answer = 100.0; glob_html_log = true; glob_good_digits = 0; glob_max_opt_iter = 10; glob_dump = false; glob_djd_debug = true; glob_display_flag = true; glob_djd_debug2 = true; glob_h_reason = 0; glob_sec_in_minute = 60; glob_min_in_hour = 60.0; glob_hours_in_day = 24.0; glob_days_in_year = 365; glob_sec_in_hour = 3600; glob_sec_in_day = 86400; glob_sec_in_year = 31536000; glob_almost_1 = 0.9990; glob_clock_sec = 0.0; glob_clock_start_sec = 0.0; glob_not_yet_finished = true; glob_initial_pass = true; glob_not_yet_start_msg = true; glob_reached_optimal_h = false; glob_optimal_done = false; glob_disp_incr = 0.1; glob_h = 0.1; glob_diff_rc_fm = 0.1; glob_diff_rc_fmm1 = 0.1; glob_diff_rc_fmm2 = 0.1; glob_diff_ord_fm = 0.1; glob_diff_ord_fmm1 = 0.1; glob_diff_ord_fmm2 = 0.1; glob_six_term_ord_save = 0.1; glob_guess_error_rc = 0.1; glob_guess_error_ord = 0.1; glob_max_h = 0.1; glob_min_h = 0.000001; glob_type_given_pole = 0; glob_large_float = 1.0e100; glob_larger_float = 1.1e100; glob_least_given_sing = 9.9e100; glob_least_ratio_sing = 9.9e100; glob_least_3_sing = 9.9e100; glob_least_6_sing = 9.9e100; glob_last_good_h = 0.1; glob_look_poles = false; glob_display_interval = 0.0; glob_next_display = 0.0; glob_dump_analytic = false; glob_abserr = 0.1e-10; glob_relerr = 0.1e-10; glob_min_pole_est = 0.1e+10; glob_max_hours = 0.0; glob_max_iter = 1000; glob_max_rel_trunc_err = 0.1e-10; glob_max_trunc_err = 0.1e-10; glob_no_eqs = 0; glob_optimal_clock_start_sec = 0.0; glob_optimal_start = 0.0; glob_upper_ratio_limit = 1.0001; glob_lower_ratio_limit = 0.9999; glob_small_float = 0.0; glob_smallish_float = 0.0; glob_unchanged_h_cnt = 0; glob_warned = false; glob_warned2 = false; glob_max_sec = 10000.0; glob_orig_start_sec = 0.0; glob_start = 0; glob_curr_iter_when_opt = 0; glob_current_iter = 0; glob_iter = 0; glob_normmax = 0.0; glob_max_minutes = 0.0; // before generate factorials init /* Initing Factorial Tables */ iiif = 0; while (iiif <= MAX_TERMS) { /* do number 1*/ jjjf = 0; while (jjjf <= MAX_TERMS) { /* do number 2*/ array_fact_1[iiif] = 0; array_fact_2[iiif][jjjf] = 0; jjjf = jjjf + 1; }/* end do number 2*/ ; iiif = iiif + 1; }/* end do number 1*/ ; /* Done Initing Factorial Table */ // before generate set diff initial array_y_set_initial[1][1] = true; array_y_set_initial[1][2] = false; array_y_set_initial[1][3] = false; array_y_set_initial[1][4] = false; array_y_set_initial[1][5] = false; array_y_set_initial[1][6] = false; array_y_set_initial[1][7] = false; array_y_set_initial[1][8] = false; array_y_set_initial[1][9] = false; array_y_set_initial[1][10] = false; array_y_set_initial[1][11] = false; array_y_set_initial[1][12] = false; array_y_set_initial[1][13] = false; array_y_set_initial[1][14] = false; array_y_set_initial[1][15] = false; array_y_set_initial[1][16] = false; array_y_set_initial[1][17] = false; array_y_set_initial[1][18] = false; array_y_set_initial[1][19] = false; array_y_set_initial[1][20] = false; array_y_set_initial[1][21] = false; array_y_set_initial[1][22] = false; array_y_set_initial[1][23] = false; array_y_set_initial[1][24] = false; array_y_set_initial[1][25] = false; array_y_set_initial[1][26] = false; array_y_set_initial[1][27] = false; array_y_set_initial[1][28] = false; array_y_set_initial[1][29] = false; array_y_set_initial[1][30] = false; array_y_set_initial[1][31] = false; array_y_set_initial[1][32] = false; array_y_set_initial[1][33] = false; array_y_set_initial[1][34] = false; array_y_set_initial[1][35] = false; array_y_set_initial[1][36] = false; array_y_set_initial[1][37] = false; array_y_set_initial[1][38] = false; array_y_set_initial[1][39] = false; array_y_set_initial[1][40] = false; // before generate init omniout const /* set default block */ /* Write Set Defaults */ glob_orig_start_sec = elapsed_time_seconds(); glob_curr_iter_when_opt = 0; glob_display_flag = true; glob_no_eqs = 1; glob_iter = -1; opt_iter = -1; glob_max_iter = 50000; glob_max_hours = 0.0; glob_max_minutes = 15.0; omniout_str(ALWAYS,"##############ECHO OF PROBLEM#################"); omniout_str(ALWAYS,"##############temp/mult_c_sinpostode.ode#################"); omniout_str(ALWAYS,"diff ( y , x , 1 ) = 2.0 * sin ( x ) ; "); omniout_str(ALWAYS,"!"); omniout_str(ALWAYS,"/* BEGIN FIRST INPUT BLOCK */"); omniout_str(ALWAYS,"Digits=32;"); omniout_str(ALWAYS,"max_terms=40;"); omniout_str(ALWAYS,"!"); omniout_str(ALWAYS,"/* END FIRST INPUT BLOCK */"); omniout_str(ALWAYS,"/* BEGIN SECOND INPUT BLOCK */"); omniout_str(ALWAYS,"x_start=0.1;"); omniout_str(ALWAYS,"x_end=5.0;"); omniout_str(ALWAYS,"array_y_init[0 + 1] = exact_soln_y(x_start);"); omniout_str(ALWAYS,"glob_look_poles=true;"); omniout_str(ALWAYS,""); omniout_str(ALWAYS,""); omniout_str(ALWAYS,""); omniout_str(ALWAYS,""); omniout_str(ALWAYS,""); omniout_str(ALWAYS,""); omniout_str(ALWAYS,""); omniout_str(ALWAYS,"glob_type_given_pole=3;"); omniout_str(ALWAYS,""); omniout_str(ALWAYS,"/* END SECOND INPUT BLOCK */"); omniout_str(ALWAYS,"/* BEGIN OVERRIDE BLOCK */"); omniout_str(ALWAYS,"glob_desired_digits_correct=10;"); omniout_str(ALWAYS,"glob_display_interval=0.01;"); omniout_str(ALWAYS,"glob_look_poles=true;"); omniout_str(ALWAYS,"glob_max_iter=100000000;"); omniout_str(ALWAYS,"glob_max_minutes=10.0;"); omniout_str(ALWAYS,"glob_subiter_method=3;"); omniout_str(ALWAYS,"/* END OVERRIDE BLOCK */"); omniout_str(ALWAYS,"!"); omniout_str(ALWAYS,"/* BEGIN USER DEF BLOCK */"); omniout_str(ALWAYS,"double exact_soln_y (double x) { "); omniout_str(ALWAYS,"return(-cos(x) * 2.0);"); omniout_str(ALWAYS,"}"); omniout_str(ALWAYS,"/* END USER DEF BLOCK */"); omniout_str(ALWAYS,"#######END OF ECHO OF PROBLEM#################"); glob_unchanged_h_cnt = 0; glob_warned = false; glob_warned2 = false; glob_small_float = 0.0; glob_smallish_float = 0.0; glob_large_float = 1.0e100; glob_larger_float = 1.1e100; glob_almost_1 = 0.99; /* before second block */ /* TOP SECOND INPUT BLOCK */ /* BEGIN SECOND INPUT BLOCK */ /* BEGIN BLOCK 2 */ /* END FIRST INPUT BLOCK */ /* BEGIN SECOND INPUT BLOCK */ x_start=0.1; x_end=5.0; array_y_init[0 + 1] = exact_soln_y(x_start); glob_look_poles=true; glob_type_given_pole=3; /* END SECOND INPUT BLOCK */ /* BEGIN OVERRIDE BLOCK */ glob_desired_digits_correct=10; glob_display_interval=0.01; glob_look_poles=true; glob_max_iter=100000000; glob_max_minutes=10.0; glob_subiter_method=3; /* END OVERRIDE BLOCK */ /* END BLOCK 2 */ /* END SECOND INPUT BLOCK */ /* BEGIN INITS AFTER SECOND INPUT BLOCK */ glob_last_good_h = glob_h; glob_max_sec = (60.0) * (glob_max_minutes) + (3600.0) * (glob_max_hours); // after second input block /* BEGIN OPTIMIZE CODE */ omniout_str(ALWAYS,"START of Optimize"); /* Start Series -- INITIALIZE FOR OPTIMIZE */ glob_check_sign = my_check_sign(x_start,x_end); found_h = false; glob_min_pole_est = glob_larger_float; last_min_pole_est = glob_larger_float; glob_least_given_sing = glob_larger_float; glob_least_ratio_sing = glob_larger_float; glob_least_3_sing = glob_larger_float; glob_least_6_sing = glob_larger_float; glob_min_h = float_abs(glob_min_h) * glob_check_sign; glob_max_h = float_abs(glob_max_h) * glob_check_sign; glob_h = float_abs(glob_min_h) * glob_check_sign; glob_display_interval = float_abs(glob_display_interval) * glob_check_sign; chk_data(); min_value = glob_larger_float; est_answer = est_size_answer(); opt_iter = 1; est_needed_step_err = estimated_needed_step_error(x_start,x_end,glob_h,est_answer); omniout_float(ALWAYS,"est_needed_step_err",32,est_needed_step_err,16,""); estimated_step_error = 0.0; while ((opt_iter <= 100) && ( ! found_h)) { /* do number 1*/ omniout_int(ALWAYS,"opt_iter",32,opt_iter,4,""); array_x[1] = x_start; array_x[2] = glob_h; glob_next_display = x_start; order_diff = 1; /* Start Series array_y */ term_no = 1; while (term_no <= order_diff) { /* do number 2*/ array_y[term_no] = array_y_init[term_no] * expt(glob_h , (term_no - 1)) / factorial_1(term_no - 1); term_no = term_no + 1; }/* end do number 2*/ ; rows = order_diff; r_order = 1; while (r_order <= rows) { /* do number 2*/ term_no = 1; while (term_no <= (rows - r_order + 1)) { /* do number 3*/ it = term_no + r_order - 1; if (term_no < MAX_TERMS) { /* if number 9*/ array_y_higher[r_order][term_no] = array_y_init[it]* expt(glob_h , (term_no - 1)) / ((factorial_1(term_no - 1))); }/* end if 9*/ ; term_no = term_no + 1; }/* end do number 3*/ ; r_order = r_order + 1; }/* end do number 2*/ ; atomall(); if (glob_check_sign * glob_max_h <= glob_check_sign * glob_h) { /* if number 9*/ omniout_str(ALWAYS,"SETTING H FOR MAX H"); glob_h = glob_check_sign * float_abs(glob_max_h); glob_h_reason = 1; found_h = true; }/* end if 9*/ ; if (glob_check_sign * glob_display_interval <= glob_check_sign * glob_h) { /* if number 9*/ omniout_str(ALWAYS,"SETTING H FOR DISPLAY INTERVAL"); glob_h_reason = 2; glob_h = glob_display_interval; found_h = true; }/* end if 9*/ ; if (glob_look_poles) { /* if number 9*/ check_for_pole(); if ((opt_iter > 2) && ( ! found_h) && ((glob_min_pole_est < 0.999 * last_min_pole_est) || (glob_min_pole_est > 1.111 * last_min_pole_est))) { /* if number 10*/ omniout_str(ALWAYS,"SETTING H FOR POLE ACCURACY"); glob_h_reason = 4; found_h = true; glob_h = glob_h/2.0; last_min_pole_est = glob_min_pole_est; } else { last_min_pole_est = glob_min_pole_est; }/* end if 10*/ ; }/* end if 9*/ ; if ( ! found_h) { /* if number 9*/ est_answer = est_size_answer(); est_needed_step_err = estimated_needed_step_error(x_start,x_end,glob_h,est_answer); omniout_float(ALWAYS,"est_needed_step_err",32,est_needed_step_err,16,""); estimated_step_error = test_suggested_h(); omniout_float(ALWAYS,"estimated_step_error",32,estimated_step_error,32,""); if (estimated_step_error < est_needed_step_err) { /* if number 10*/ omniout_str(ALWAYS,"Double H and LOOP"); glob_h = glob_h*2.0; } else { omniout_str(ALWAYS,"Found H for OPTIMAL"); found_h = true; glob_h_reason = 3; glob_h = glob_h/2.0; }/* end if 10*/ ; }/* end if 9*/ ; opt_iter = opt_iter + 1; }/* end do number 1*/ ; if (( ! found_h) && (opt_iter == 1)) { /* if number 9*/ omniout_str(ALWAYS,"Beginning glob_h too large."); found_h = false; }/* end if 9*/ ; /* END OPTIMIZE CODE */ if (glob_html_log) { /* if number 9*/ html_log_file = fopen("entry.html","w"); }/* end if 9*/ ; /* BEGIN SOLUTION CODE */ if (found_h) { /* if number 9*/ omniout_str(ALWAYS,"START of Soultion"); /* Start Series -- INITIALIZE FOR SOLUTION */ array_x[1] = x_start; array_x[2] = glob_h; glob_next_display = x_start; glob_min_pole_est = glob_larger_float; glob_least_given_sing = glob_larger_float; glob_least_ratio_sing = glob_larger_float; glob_least_3_sing = glob_larger_float; glob_least_6_sing = glob_larger_float; order_diff = 1; /* Start Series array_y */ term_no = 1; while (term_no <= order_diff) { /* do number 1*/ array_y[term_no] = array_y_init[term_no] * expt(glob_h , (term_no - 1)) / factorial_1(term_no - 1); term_no = term_no + 1; }/* end do number 1*/ ; rows = order_diff; r_order = 1; while (r_order <= rows) { /* do number 1*/ term_no = 1; while (term_no <= (rows - r_order + 1)) { /* do number 2*/ it = term_no + r_order - 1; if (term_no < MAX_TERMS) { /* if number 10*/ array_y_higher[r_order][term_no] = array_y_init[it]* expt(glob_h , (term_no - 1)) / ((factorial_1(term_no - 1))); }/* end if 10*/ ; term_no = term_no + 1; }/* end do number 2*/ ; r_order = r_order + 1; }/* end do number 1*/ ; current_iter = 1; glob_clock_start_sec = elapsed_time_seconds(); glob_clock_sec = elapsed_time_seconds(); glob_current_iter = 0; glob_iter = 0; omniout_str(DEBUGL," "); glob_reached_optimal_h = true; glob_optimal_clock_start_sec = elapsed_time_seconds(); while ((glob_current_iter < glob_max_iter) && (glob_check_sign * array_x[1] < glob_check_sign * x_end ) && (((glob_clock_sec) - (glob_orig_start_sec)) < (glob_max_sec))) { /* do number 1*/ /* left paren 0001C */ if (reached_interval()) { /* if number 10*/ omniout_str(INFO," "); omniout_str(INFO,"TOP MAIN SOLVE Loop"); }/* end if 10*/ ; glob_iter = glob_iter + 1; glob_clock_sec = elapsed_time_seconds(); track_estimated_error(); glob_current_iter = glob_current_iter + 1; atomall(); track_estimated_error(); display_alot(current_iter); check_for_pole(); if (reached_interval()) { /* if number 10*/ glob_next_display = glob_next_display + glob_display_interval; }/* end if 10*/ ; array_x[1] = array_x[1] + glob_h; array_x[2] = glob_h; /* Jump Series array_y */; order_diff = 2; /* START PART 1 SUM AND ADJUST */ /* START SUM AND ADJUST EQ =1 */ /* sum_and_adjust array_y */ /* BEFORE ADJUST SUBSERIES EQ =1 */ ord = 2; calc_term = 1; /* adjust_subseriesarray_y */ iii = MAX_TERMS; while (iii >= calc_term) { /* do number 2*/ array_y_higher_work[2][iii] = array_y_higher[2][iii] / expt(glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1); iii = iii - 1; }/* end do number 2*/ ; /* AFTER ADJUST SUBSERIES EQ =1 */ /* BEFORE SUM SUBSERIES EQ =1 */ temp_sum = 0.0; ord = 2; calc_term = 1; /* sum_subseriesarray_y */ iii = MAX_TERMS; while (iii >= calc_term) { /* do number 2*/ temp_sum = temp_sum + array_y_higher_work[ord][iii]; iii = iii - 1; }/* end do number 2*/ ; array_y_higher_work2[ord][calc_term] = temp_sum * expt(glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1)); /* AFTER SUM SUBSERIES EQ =1 */ /* BEFORE ADJUST SUBSERIES EQ =1 */ ord = 1; calc_term = 2; /* adjust_subseriesarray_y */ iii = MAX_TERMS; while (iii >= calc_term) { /* do number 2*/ array_y_higher_work[1][iii] = array_y_higher[1][iii] / expt(glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1); iii = iii - 1; }/* end do number 2*/ ; /* AFTER ADJUST SUBSERIES EQ =1 */ /* BEFORE SUM SUBSERIES EQ =1 */ temp_sum = 0.0; ord = 1; calc_term = 2; /* sum_subseriesarray_y */ iii = MAX_TERMS; while (iii >= calc_term) { /* do number 2*/ temp_sum = temp_sum + array_y_higher_work[ord][iii]; iii = iii - 1; }/* end do number 2*/ ; array_y_higher_work2[ord][calc_term] = temp_sum * expt(glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1)); /* AFTER SUM SUBSERIES EQ =1 */ /* BEFORE ADJUST SUBSERIES EQ =1 */ ord = 1; calc_term = 1; /* adjust_subseriesarray_y */ iii = MAX_TERMS; while (iii >= calc_term) { /* do number 2*/ array_y_higher_work[1][iii] = array_y_higher[1][iii] / expt(glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1); iii = iii - 1; }/* end do number 2*/ ; /* AFTER ADJUST SUBSERIES EQ =1 */ /* BEFORE SUM SUBSERIES EQ =1 */ temp_sum = 0.0; ord = 1; calc_term = 1; /* sum_subseriesarray_y */ iii = MAX_TERMS; while (iii >= calc_term) { /* do number 2*/ temp_sum = temp_sum + array_y_higher_work[ord][iii]; iii = iii - 1; }/* end do number 2*/ ; array_y_higher_work2[ord][calc_term] = temp_sum * expt(glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1)); /* AFTER SUM SUBSERIES EQ =1 */ /* END SUM AND ADJUST EQ =1 */ /* END PART 1 */ /* START PART 2 MOVE TERMS to REGULAR Array */ term_no = MAX_TERMS; while (term_no >= 1) { /* do number 2*/ array_y[term_no] = array_y_higher_work2[1][term_no]; ord = 1; while (ord <= order_diff) { /* do number 3*/ array_y_higher[ord][term_no] = array_y_higher_work2[ord][term_no]; ord = ord + 1; }/* end do number 3*/ ; term_no = term_no - 1; }/* end do number 2*/ ; /* END PART 2 HEVE MOVED TERMS to REGULAR Array */ ; }/* end do number 1*/ ;/* right paren 0001C */ omniout_str(ALWAYS,"Finished!"); if (glob_iter >= glob_max_iter) { /* if number 10*/ omniout_str(ALWAYS,"Maximum Iterations Reached before Solution Completed!"); }/* end if 10*/ ; if (elapsed_time_seconds() - (glob_orig_start_sec) >= (glob_max_sec )) { /* if number 10*/ omniout_str(ALWAYS,"Maximum Time Reached before Solution Completed!"); }/* end if 10*/ ; glob_clock_sec = elapsed_time_seconds(); omniout_str(INFO,"diff ( y , x , 1 ) = 2.0 * sin ( x ) ; "); omniout_int(INFO,"Iterations ",32,glob_iter,4," ") ; prog_report(x_start,x_end); if (glob_html_log) { /* if number 10*/ logstart(html_log_file); logitem_str(html_log_file,"2014-01-30T02:06:01-06:00") ; logitem_str(html_log_file,"c") ; logitem_str(html_log_file,"mult_c_sin") ; logitem_str(html_log_file,"diff ( y , x , 1 ) = 2.0 * sin ( x ) ; ") ; logitem_float(html_log_file,x_start) ; logitem_float(html_log_file,x_end) ; logitem_float(html_log_file,array_x[1]) ; logitem_float(html_log_file,glob_h) ; logitem_h_reason(html_log_file) ; logitem_str(html_log_file,"16") ; logitem_float(html_log_file,glob_desired_digits_correct) ; if (array_est_digits[1] != -16) { /* if number 11*/ logitem_integer(html_log_file,array_est_digits[1]) ; } else { logitem_str(html_log_file,"Unknown") ; }/* end if 11*/ if (array_good_digits[1] != -16) { /* if number 11*/ logitem_integer(html_log_file,array_good_digits[1]) ; } else { logitem_str(html_log_file,"Unknown") ; }/* end if 11*/ logitem_integer(html_log_file,MAX_TERMS) ; if (glob_type_given_pole == 0) { /* if number 11*/ logitem_str(html_log_file,"Not Given") ; logitem_str(html_log_file,"NA") ; } else if (glob_type_given_pole == 4) { /* if number 12*/ logitem_str(html_log_file,"No Solution") ; logitem_str(html_log_file,"NA") ; } else if (glob_type_given_pole == 5) { /* if number 13*/ logitem_str(html_log_file,"Some Pole") ; logitem_str(html_log_file,"????") ; } else if (glob_type_given_pole == 3) { /* if number 14*/ logitem_str(html_log_file,"No Pole") ; logitem_str(html_log_file,"NA") ; } else if (glob_type_given_pole == 1) { /* if number 15*/ logitem_str(html_log_file,"Real Sing") ; logitem_float(html_log_file,glob_least_given_sing) ; } else if (glob_type_given_pole == 2) { /* if number 16*/ logitem_str(html_log_file,"Complex Sing") ; logitem_float(html_log_file,glob_least_given_sing) ; }/* end if 16*/ ; if (glob_least_ratio_sing < glob_large_float) { /* if number 16*/ logitem_float(html_log_file,glob_least_ratio_sing) ; } else { logitem_str(html_log_file,"NONE") ; }/* end if 16*/ ; if (glob_least_3_sing < glob_large_float) { /* if number 16*/ logitem_float(html_log_file,glob_least_3_sing) ; } else { logitem_str(html_log_file,"NONE") ; }/* end if 16*/ ; if (glob_least_6_sing < glob_large_float) { /* if number 16*/ logitem_float(html_log_file,glob_least_6_sing) ; } else { logitem_str(html_log_file,"NONE") ; }/* end if 16*/ ; logitem_integer(html_log_file,glob_iter) ; logitem_time(html_log_file,(glob_clock_sec)) ; if (glob_percent_done < 100.0) { /* if number 16*/ logitem_time(html_log_file,(glob_total_exp_sec)) ; 0; } else { logitem_str(html_log_file,"Done") ; 0; }/* end if 16*/ ; log_revs(html_log_file," 255 ") ; logitem_str(html_log_file,"mult_c_sin diffeq.c") ; logitem_str(html_log_file,"mult_c_sin c results") ; logitem_str(html_log_file,"OK") ; logend(html_log_file) ; ; }/* end if 15*/ ; if (glob_html_log) { /* if number 15*/ fclose(html_log_file); }/* end if 15*/ ; ;; }/* end if 14*/ /* END OUTFILEMAIN */ } /* End Function number1 */