/* BEGIN OUTFILE1 */ #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 30 #include #include #include #include #define convfp(x) ((double)x) #define convfloat(x) ((double)x) int glob_max_terms; long elapsed_time_seconds(); int reached_interval(); int not_reached_end(double ,double ); double expt(double ,double ); double arcsin(double ); double arccos(double ); double arctan(double ); double omniabs(double ); double ln(double ); double Si(double ); double Ci(double ); double ats(int,double *,double *,int); double att(int,double *,double *,int); double factorial_1(int); double factorial_3(int,int); void display_pole(); 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,int secs_in); void omniout_timestr(int secs_in); double array_const_m1[max_terms + 1]; 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 check_sign(double x0,double xf); double exact_soln_y(double ); /* Top Generate Globals Definition */ int MAX_UNCHANGED=10; double glob_check_sign=1.0; double glob_desired_digits_correct=8.0; double glob_max_value3=0.0; double glob_ratio_of_radius=0.01; 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; 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_sec_in_minute=60; int glob_min_in_hour=60; int glob_hours_in_day=24; 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_max_h=0.1; double glob_large_float=9.0e100; double glob_last_good_h=0.1; unsigned char glob_look_poles=false; unsigned char glob_neg_h=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_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_small_float=0.1e-200; double glob_smallish_float=0.1e-100; 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 */ double array_y_init[max_terms + 1]; double array_norms[max_terms + 1]; double array_fact_1[max_terms + 1]; double array_pole[max_terms + 1]; double array_1st_rel_error[max_terms + 1]; double array_last_rel_error[max_terms + 1]; double array_type_pole[max_terms + 1]; double array_y[max_terms + 1]; double array_x[max_terms + 1]; double array_tmp0[max_terms + 1]; double array_tmp1[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_poles[1 + 1][3+ 1]; double array_real_pole[1 + 1][3+ 1]; double array_complex_pole[1 + 1][3+ 1]; double array_fact_2[max_terms + 1][max_terms+ 1]; double array_const_1[max_terms+1]; double array_const_0D0[max_terms+1]; double check_sign(double x0 ,double xf) { double ret; if (xf > x0) { /* if number 1*/ ret = 1.0; } else { ret = -1.0; }/* end if 1*/ ; return ret;; } /* End Function number1 */ double est_size_answer() { double min_size; min_size = glob_large_float; if (omniabs(array_y[1]) < min_size) { /* if number 1*/ min_size = omniabs(array_y[1]); omniout_float(ALWAYS,"min_size",32,min_size,32,""); }/* end if 1*/ ; if (min_size < 1.0) { /* if number 1*/ min_size = 1.0; omniout_float(ALWAYS,"min_size",32,min_size,32,""); }/* end if 1*/ ; return min_size; } /* End Function number1 */ double test_suggested_h() { double max_value3,hn_div_ho,hn_div_ho_2,hn_div_ho_3,value3; int no_terms; max_value3 = 0.0; no_terms = glob_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,""); value3 = omniabs(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 (value3 > max_value3) { /* if number 1*/ max_value3 = value3; omniout_float(ALWAYS,"value3",32,value3,32,""); }/* end if 1*/ ; omniout_float(ALWAYS,"max_value3",32,max_value3,32,""); return max_value3; } /* End Function number1 */ int reached_interval() { int ret; if (glob_check_sign * (array_x[1]) >= glob_check_sign * glob_next_display) { /* if number 1*/ ret = true; } else { ret = false; }/* end if 1*/ ; return(ret); } /* End Function number1 */ void display_alot(int iter) { double abserr, analytic_val_y, ind_var, numeric_val, relerr; int term_no; /* TOP DISPLAY ALOT */ if (reached_interval()) { /* if number 1*/ if (iter >= 0) { /* if number 2*/ 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 = omniabs(numeric_val - analytic_val_y); omniout_float(ALWAYS,"y[1] (numeric) ",33,numeric_val,20," "); if (omniabs(analytic_val_y) != 0.0) { /* if number 3*/ relerr = abserr*100.0/omniabs(analytic_val_y); if (relerr > 0.0000000000000000000000000000000001) { /* if number 4*/ glob_good_digits = -trunc(log10(relerr)) + 2; } else { glob_good_digits = 16; }/* end if 4*/ ; } else { relerr = -1.0 ; glob_good_digits = -1; }/* end if 3*/ ; if (glob_iter == 1) { /* if number 3*/ array_1st_rel_error[1] = relerr; } else { array_last_rel_error[1] = relerr; }/* end if 3*/ ; omniout_float(ALWAYS,"absolute error ",4,abserr,20," "); omniout_float(ALWAYS,"relative error ",4,relerr,20,"%"); omniout_int(INFO,"Correct digits ",32,glob_good_digits,4," ") ; omniout_float(ALWAYS,"h ",4,glob_h,20," "); }/* end if 2*/ ; /* BOTTOM DISPLAY ALOT */ }/* end if 1*/ ; } /* End Function number1 */ double adjust_for_pole(double h_param) { double hnew, sz2, tmp; /* TOP ADJUST FOR POLE */ hnew = h_param; glob_normmax = glob_small_float; if (omniabs(array_y_higher[1][1]) > glob_small_float) { /* if number 1*/ tmp = omniabs(array_y_higher[1][1]); if (tmp < glob_normmax) { /* if number 2*/ glob_normmax = tmp; }/* end if 2*/ }/* end if 1*/ ; if (glob_look_poles && (omniabs(array_pole[1]) > glob_small_float) && (array_pole[1] != glob_large_float)) { /* if number 1*/ sz2 = array_pole[1]/10.0; if (sz2 < hnew) { /* if number 2*/ omniout_float(INFO,"glob_h adjusted to ",20,h_param,12,"due to singularity."); omniout_str(INFO,"Reached Optimal"); return(hnew); }/* end if 2*/ }/* end if 1*/ ; if ( ! glob_reached_optimal_h) { /* 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]; }/* end if 1*/ ; hnew = sz2; ;/* END block */ return(hnew); /* BOTTOM ADJUST FOR POLE */ } /* End Function number1 */ void prog_report(double x_start,double x_end) { int 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 = 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)); glob_total_exp_sec = glob_optimal_expect_sec + total_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)) { /* 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)); omniout_str_noeol(INFO,"Expected Total Time "); omniout_timestr(convfloat(glob_total_exp_sec)); }/* 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 number1 */ void check_for_pole() { int cnt, m, n, found_sing, term; double dr1, dr2, ds1, ds2, hdrc,hdrc_BBB, nr1, nr2, ord_no, rad_c, rcs, rm0, rm1, rm2, rm3, rm4, h_new, ratio; /* 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) && ((omniabs(array_y_higher[1][m]) < glob_small_float * glob_small_float) || (omniabs(array_y_higher[1][m-1]) < glob_small_float * glob_small_float) || (omniabs(array_y_higher[1][m-2]) < glob_small_float * glob_small_float ))) { /* do number 2*/ m = m - 1; }/* end do number 2*/ ; if (m > 10) { /* 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)*rm0-convfloat(m-1)*rm1; if (omniabs(hdrc) > glob_small_float * glob_small_float) { /* if number 2*/ rcs = glob_h/hdrc; ord_no = (rm1*convfloat((m-2)*(m-2))-rm0*convfloat(m-3))/hdrc; 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; }/* end if 2*/ } else { array_real_pole[1][1] = glob_large_float; array_real_pole[1][2] = glob_large_float; }/* 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) && (n >= 10)) { /* do number 2*/ if (omniabs(array_y_higher[1][n]) > glob_small_float) { /* if number 1*/ cnt = cnt + 1; } else { cnt = 0; }/* end if 1*/ ; n = n - 1; }/* end do number 2*/ ; m = n + cnt; if (m <= 10) { /* if number 1*/ rad_c = glob_large_float; ord_no = glob_large_float; } else if (((omniabs(array_y_higher[1][m]) >= (glob_large_float)) || (omniabs(array_y_higher[1][m-1]) >=(glob_large_float)) || (omniabs(array_y_higher[1][m-2]) >= (glob_large_float)) || (omniabs(array_y_higher[1][m-3]) >= (glob_large_float)) || (omniabs(array_y_higher[1][m-4]) >= (glob_large_float)) || (omniabs(array_y_higher[1][m-5]) >= (glob_large_float))) || ((omniabs(array_y_higher[1][m]) <= (glob_small_float)) || (omniabs(array_y_higher[1][m-1]) <=(glob_small_float)) || (omniabs(array_y_higher[1][m-2]) <= (glob_small_float)) || (omniabs(array_y_higher[1][m-3]) <= (glob_small_float)) || (omniabs(array_y_higher[1][m-4]) <= (glob_small_float)) || (omniabs(array_y_higher[1][m-5]) <= (glob_small_float)))) { /* if number 2*/ rad_c = glob_large_float; ord_no = 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 ((omniabs(nr1 * dr2 - nr2 * dr1) <= glob_small_float) || (omniabs(dr1) <= glob_small_float)) { /* if number 3*/ rad_c = glob_large_float; ord_no = glob_large_float; } else { if (omniabs(nr1*dr2 - nr2 * dr1) > glob_small_float) { /* 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 (omniabs(rcs) > glob_small_float) { /* if number 5*/ if (rcs > 0.0) { /* if number 6*/ rad_c = sqrt(rcs) * omniabs(glob_h); } else { rad_c = glob_large_float; }/* end if 6*/ } else { rad_c = glob_large_float; ord_no = glob_large_float; }/* end if 5*/ } else { rad_c = glob_large_float; ord_no = glob_large_float; }/* end if 4*/ }/* end if 3*/ ; array_complex_pole[1][1] = rad_c; array_complex_pole[1][2] = ord_no; }/* end if 2*/ ; /* BOTTOM RADII COMPLEX EQ = 1 */ found_sing = 0; /* TOP WHICH RADII EQ = 1 */ if (1 != found_sing && ((array_real_pole[1][1] == glob_large_float) || (array_real_pole[1][2] == glob_large_float)) && ((array_complex_pole[1][1] != glob_large_float) && (array_complex_pole[1][2] != glob_large_float)) && ((array_complex_pole[1][1] > 0.0) && (array_complex_pole[1][2] > 0.0))) { /* if number 2*/ array_poles[1][1] = array_complex_pole[1][1]; array_poles[1][2] = array_complex_pole[1][2]; found_sing = 1; array_type_pole[1] = 2; if (glob_display_flag) { /* if number 3*/ if (reached_interval()) { /* if number 4*/ omniout_str(ALWAYS,"Complex estimate of poles used for equation 1"); }/* end if 4*/ ; }/* end if 3*/ ; }/* end if 2*/ ; if (1 != found_sing && ((array_real_pole[1][1] != glob_large_float) && (array_real_pole[1][2] != glob_large_float) && (array_real_pole[1][1] > 0.0) && (array_real_pole[1][2] > -1.0 * glob_smallish_float) && ((array_complex_pole[1][1] == glob_large_float) || (array_complex_pole[1][2] == glob_large_float) || (array_complex_pole[1][1] <= 0.0 ) || (array_complex_pole[1][2] <= 0.0)))) { /* if number 2*/ array_poles[1][1] = array_real_pole[1][1]; array_poles[1][2] = array_real_pole[1][2]; found_sing = 1; array_type_pole[1] = 1; if (glob_display_flag) { /* if number 3*/ if (reached_interval()) { /* if number 4*/ omniout_str(ALWAYS,"Real estimate of pole used for equation 1"); }/* end if 4*/ ; }/* end if 3*/ ; }/* end if 2*/ ; if (1 != found_sing && (((array_real_pole[1][1] == glob_large_float) || (array_real_pole[1][2] == glob_large_float)) && ((array_complex_pole[1][1] == glob_large_float) || (array_complex_pole[1][2] == glob_large_float)))) { /* if number 2*/ array_poles[1][1] = glob_large_float; array_poles[1][2] = glob_large_float; found_sing = 1; array_type_pole[1] = 3; if (reached_interval()) { /* if number 3*/ omniout_str(ALWAYS,"NO POLE for equation 1"); }/* end if 3*/ ; }/* end if 2*/ ; if (1 != found_sing && ((array_real_pole[1][1] < array_complex_pole[1][1]) && (array_real_pole[1][1] > 0.0) && (array_real_pole[1][2] > -1.0 * glob_smallish_float))) { /* if number 2*/ array_poles[1][1] = array_real_pole[1][1]; array_poles[1][2] = array_real_pole[1][2]; found_sing = 1; array_type_pole[1] = 1; if (glob_display_flag) { /* if number 3*/ if (reached_interval()) { /* if number 4*/ omniout_str(ALWAYS,"Real estimate of pole used for equation 1"); }/* end if 4*/ ; }/* end if 3*/ ; }/* end if 2*/ ; if (1 != found_sing && ((array_complex_pole[1][1] != glob_large_float) && (array_complex_pole[1][2] != glob_large_float) && (array_complex_pole[1][1] > 0.0) && (array_complex_pole[1][2] > 0.0))) { /* 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_sing = 1; if (glob_display_flag) { /* if number 3*/ if (reached_interval()) { /* if number 4*/ omniout_str(ALWAYS,"Complex estimate of poles used for equation 1"); }/* end if 4*/ ; }/* end if 3*/ ; }/* end if 2*/ ; if (1 != found_sing ) { /* if number 2*/ array_poles[1][1] = glob_large_float; array_poles[1][2] = glob_large_float; array_type_pole[1] = 3; if (reached_interval()) { /* if number 3*/ omniout_str(ALWAYS,"NO POLE for equation 1"); }/* end if 3*/ ; }/* 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]) { /* if number 2*/ array_pole[1] = array_poles[1][1]; array_pole[2] = array_poles[1][2]; }/* end if 2*/ ; /* BOTTOM WHICH RADIUS EQ = 1 */ /* START ADJUST ALL SERIES */ if (array_pole[1] * glob_ratio_of_radius < omniabs(glob_h)) { /* if number 2*/ h_new = array_pole[1] * glob_ratio_of_radius; term = 1; ratio = 1.0; while (term <= glob_max_terms) { /* do number 2*/ 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 / omniabs(glob_h); term = term + 1; }/* end do number 2*/ ; glob_h = h_new; }/* end if 2*/ ; /* BOTTOM ADJUST ALL SERIES */ if (reached_interval()) { /* if number 2*/ display_pole(); }/* end if 2*/ } /* End Function number1 */ void get_norms() { int iii; if ( ! glob_initial_pass) { /* if number 2*/ iii = 1; while (iii <= glob_max_terms) { /* do number 2*/ array_norms[iii] = 0.0; iii = iii + 1; }/* end do number 2*/ ; /* TOP GET NORMS */ iii = 1; while (iii <= glob_max_terms) { /* do number 2*/ if (omniabs(array_y[iii]) > array_norms[iii]) { /* if number 3*/ array_norms[iii] = omniabs(array_y[iii]); }/* end if 3*/ ; iii = iii + 1; }/* end do number 2*/ /* BOTTOM GET NORMS */ ; }/* end if 2*/ ; } /* End Function number1 */ void atomall() { int kkk, order_d, term, adj2, adj3; double temporary,temp,temp2; /* TOP ATOMALL */ /* END OUTFILE1 */ /* BEGIN ATOMHDR1 */ /* emit pre add CONST FULL $eq_no = 1 i = 1 */ array_tmp1[1] = array_const_0D0[1] + array_y[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 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[1] * expt(glob_h , (1)) * factorial_3(0,1); array_y[2] = temporary; array_y_higher[1][2] = temporary; 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 add CONST FULL $eq_no = 1 i = 2 */ array_tmp1[2] = array_y[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 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[2] * expt(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; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 3; /* END ATOMHDR2 */ /* BEGIN ATOMHDR3 */ /* emit pre add CONST FULL $eq_no = 1 i = 3 */ array_tmp1[3] = array_y[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 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[3] * expt(glob_h , (1)) * factorial_3(2,3); array_y[4] = temporary; array_y_higher[1][4] = temporary; 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 add CONST FULL $eq_no = 1 i = 4 */ array_tmp1[4] = array_y[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 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[4] * expt(glob_h , (1)) * factorial_3(3,4); array_y[5] = temporary; array_y_higher[1][5] = temporary; 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 add CONST FULL $eq_no = 1 i = 5 */ array_tmp1[5] = array_y[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 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[5] * expt(glob_h , (1)) * factorial_3(4,5); array_y[6] = temporary; array_y_higher[1][6] = temporary; 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 <= glob_max_terms) { /* do number 1*/ /* END OUTFILE3 */ /* BEGIN OUTFILE4 */ /* emit NOT FULL - FULL add $eq_no = 1 */ array_tmp1[kkk] = array_y[kkk]; /* emit assign $eq_no = 1 */ order_d = 1; if (kkk + order_d + 1 <= glob_max_terms) { /* if number 1*/ if ( ! array_y_set_initial[1][kkk + order_d]) { /* if number 2*/ temporary = array_tmp1[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) { /* do number 2*/ if (adj3 <= order_d + 1) { /* if number 3*/ if (adj2 > 0) { /* if number 4*/ temporary = temporary / glob_h * convfp(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 2*/ }/* end if 2*/ }/* end if 1*/ ; kkk = kkk + 1; }/* end do number 1*/ ; /* BOTTOM ATOMALL */ /* END OUTFILE4 */ /* BEGIN OUTFILE5 */ /* BOTTOM ATOMALL ??? */ } /* End Function number1 */ /* 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 dump_series(int iolevel,char *dump_label,char *series_name,double *arr_series,int numb) { int i; if (glob_iolevel >= iolevel) { i = 1; while (i <= numb) { printf("%s %s [ %d ] %g\n" , dump_label,series_name ,i,arr_series[i]); i++; } } } /* End Function number1 */ void cs_info(int iolevel,char *str) { if (glob_iolevel >= iolevel) { /* if number 1*/ printf("cs_info %s glob_h: = %g\n", str,glob_h); }/* end if 1*/ } /* End Function number1 */ void logitem_time(FILE *fd,int secs_in) { int days_int, hours_int, minutes_int, sec_int, years_int,sec_temp; fprintf(fd,""); if (secs_in >= 0) { /* if number 1*/ years_int = (secs_in / glob_sec_in_year); sec_temp = (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); if (years_int > 0) { /* if number 2*/ fprintf(fd,"%d Years %d Days %d Hours %d Minutes %d Seconds",years_int,days_int,hours_int,minutes_int,sec_int); } else if (days_int > 0) { /* if number 3*/ fprintf(fd,"%d Days %d Hours %d Minutes %d Seconds",days_int,hours_int,minutes_int,sec_int); } else if (hours_int > 0) { /* if number 4*/ fprintf(fd,"%d Hours %d Minutes %d Seconds",hours_int,minutes_int,sec_int); } else if (minutes_int > 0) { /* if number 5*/ fprintf(fd,"%d Minutes %d Seconds",minutes_int,sec_int); } else { fprintf(fd,"%d Seconds",sec_int); }/* end if 5*/ } else { fprintf(fd,"Unknown"); }/* end if 4*/ fprintf(fd,""); } /* End Function number1 */ void omniout_timestr(int secs_in) { int days_int, hours_int, minutes_int, sec_int, years_int,sec_temp; if (secs_in >= 0) { /* if number 4*/ years_int = (secs_in / glob_sec_in_year); sec_temp = (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); if (years_int > 0) { /* if number 5*/ printf(" = %d Years %d Days %d Hours %d Minutes %d Seconds\n",years_int,days_int,hours_int,minutes_int,sec_int); } else if (days_int > 0) { /* if number 6*/ printf(" = %d Days %d Hours %d Minutes %d Seconds\n",days_int,hours_int,minutes_int,sec_int); } else if (hours_int > 0) { /* if number 7*/ printf(" = %d Hours %d Minutes %d Seconds\n",hours_int,minutes_int,sec_int); } else if (minutes_int > 0) { /* if number 8*/ printf(" = %d Minutes %d Seconds\n",minutes_int,sec_int); } else { printf(" = %d Seconds\n",sec_int); }/* end if 8*/ } else { printf(" Unknown\n"); }/* end if 7*/ } /* End Function number1 */ double ats(int mmm_ats,double *arr_a,double *arr_b,int jjj_ats) { int iii_ats, lll_ats; double ma_ats, 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; ret_ats = ret_ats + arr_a[iii_ats]*arr_b[lll_ats]; 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) { /* do number 1*/ lll_att = ma_att - iii_att; al_att = (lll_att - 1); if (lll_att <= glob_max_terms) { /* if number 8*/ ret_att = ret_att + arr_aa[iii_att]*arr_bb[lll_att]* convfp(al_att); }/* end if 8*/ ; iii_att = iii_att + 1; }/* end do number 1*/ ; ret_att = ret_att / convfp(mmm_att) ; }/* end if 7*/ ; return ret_att; } /* End Function number1 */ void display_pole_debug(int typ,double radius,double order2) { if (typ == 1) { /* if number 7*/ omniout_str(ALWAYS,"Real"); } else { omniout_str(ALWAYS,"Complex"); }/* end if 7*/ ; omniout_float(ALWAYS,"DBG Radius of convergence ",4, radius,4," "); omniout_float(ALWAYS,"DBG Order of pole ",4, order2,4," "); } /* End Function number1 */ void display_pole() { if ((array_pole[1] != glob_large_float) && (array_pole[1] > 0.0) && (array_pole[2] != glob_large_float) && (array_pole[2]> 0.0) && glob_display_flag) { /* if number 7*/ omniout_float(ALWAYS,"Radius of convergence ",4, array_pole[1],4," "); omniout_float(ALWAYS,"Order of pole ",4, array_pole[2],4," "); }/* end if 7*/ } /* 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.0000000000000000000000000000000001) { /* if number 8*/ good_digits = 1-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_pole(FILE * file,int pole) { fprintf(file,""); if (pole == 0) { /* if number 7*/ fprintf(file,"NA"); } else if (pole == 1) { /* if number 8*/ fprintf(file,"Real"); } else if (pole == 2) { /* if number 9*/ fprintf(file,"Complex"); } else { fprintf(file,"No Pole"); }/* end if 9*/ 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_terms < 15) || (glob_max_terms > 512)) { /* if number 9*/ omniout_str(ALWAYS,"Illegal max_terms = -- Using 30"); glob_max_terms = 30; }/* end if 9*/ ; if (glob_max_iter < 2) { /* if number 9*/ omniout_str(ALWAYS,"Illegal max_iter"); errflag = true; }/* end if 9*/ ; if (errflag) { /* if number 9*/ }/* end if 9*/ } /* 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 9*/ sec_left = 0.0; } else { if (sub2 > 0.0) { /* if number 10*/ rrr = (sub1/sub2); sec_left = rrr * ms2 - ms2; } else { sec_left = 0.0; }/* end if 10*/ }/* end if 9*/ ; 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 9*/ rrr = (100.0*sub2)/sub1; } else { rrr = 0.0; }/* end if 9*/ ; return rrr; } /* End Function number1 */ double factorial_2(int nnn) { double ret; if (nnn > 0) { /* if number 9*/ return nnn * factorial_2(nnn - 1); } else { return 1.0; }/* end if 9*/ } /* End Function number1 */ double factorial_1(int nnn) { double ret; if (nnn <= glob_max_terms) { /* if number 9*/ if (array_fact_1[nnn] == 0) { /* if number 10*/ ret = factorial_2(nnn); array_fact_1[nnn] = ret; } else { ret = array_fact_1[nnn]; }/* end if 10*/ ; } else { ret = factorial_2(nnn); }/* end if 9*/ ; return ret; } /* End Function number1 */ double factorial_3(int mmm,int nnn) { double ret; if ((nnn <= glob_max_terms) && (mmm <= glob_max_terms)) { /* if number 9*/ if (array_fact_2[mmm][nnn] == 0) { /* if number 10*/ ret = factorial_1(mmm)/factorial_1(nnn); array_fact_2[mmm][nnn] = ret; } else { ret = array_fact_2[mmm][nnn]; }/* end if 10*/ ; } else { ret = factorial_2(mmm)/factorial_2(nnn); }/* end if 9*/ ; return ret; } /* End Function number1 */ long elapsed_time_seconds() { struct timeval t; struct timezone tz; gettimeofday(&t,&tz); return (t.tv_sec); } /* End Function number1 */ double expt(double x,double y) { return pow(x,y); } /* End Function number1 */ double arcsin(double x) { return asin(x); } /* End Function number1 */ double arccos(double x) { return acos(x); } /* End Function number1 */ double arctan(double x) { return atan(x); } /* End Function number1 */ double omniabs(double x) { return fabs(x); } /* End Function number1 */ double ln(double x) { return log(x); } /* End Function number1 */ double Si(double x) { return (0.0); } /* End Function number1 */ double Ci(double x) { return (0.0); } /* 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) * omniabs(estimated_answer); 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 = omniabs(desired_abs_gbl_error / estimated_steps); omniout_float(ALWAYS,"step_error",32,step_error,32,""); return (step_error);; } /* End Function number1 */ /* END ATS LIBRARY BLOCK */ /* BEGIN USER DEF BLOCK */ /* BEGIN USER DEF BLOCK */ double exact_soln_y (double x) { return(exp(x)); } /* END USER DEF BLOCK */ /* END USER DEF BLOCK */ /* END OUTFILE5 */ int main() { /* BEGIN OUTFIEMAIN */ double d1,d2,d3,d4,est_err_2,Digits; int niii,done_once,it,opt_iter; 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; double temp_sum, est_needed_step_err,value3,min_value,est_answer,best_h,found_h; double x_start;double x_end; double tmp; /* Write Set Defaults */ glob_orig_start_sec = elapsed_time_seconds(); MAX_UNCHANGED = 10; 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/diff0postode.ode#################"); omniout_str(ALWAYS,"diff ( y , x , 1 ) = y;"); omniout_str(ALWAYS,"!"); omniout_str(ALWAYS,"/* BEGIN FIRST INPUT BLOCK */"); omniout_str(ALWAYS,"Digits=32;"); omniout_str(ALWAYS,"max_terms=30;"); omniout_str(ALWAYS,"!"); omniout_str(ALWAYS,"/* END FIRST INPUT BLOCK */"); omniout_str(ALWAYS,"/* BEGIN SECOND INPUT BLOCK */"); omniout_str(ALWAYS,"x_start=-5.0;"); 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=false;"); omniout_str(ALWAYS,"glob_max_iter=1000000;"); omniout_str(ALWAYS,"glob_display_interval=0.1;"); omniout_str(ALWAYS,"glob_max_minutes=10;"); 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.001;"); omniout_str(ALWAYS,"glob_look_poles=true;"); omniout_str(ALWAYS,"glob_max_iter=10000000;"); omniout_str(ALWAYS,"glob_max_minutes=3;"); 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(exp(x));"); omniout_str(ALWAYS,"}"); 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 = 1.0e-200; glob_smallish_float = 1.0e-64; glob_large_float = 1.0e100; glob_almost_1 = 0.99; /* BEGIN FIRST INPUT BLOCK */ /* END FIRST INPUT BLOCK */ /* START OF INITS AFTER INPUT BLOCK */ glob_max_terms = max_terms; glob_html_log = true; /* END OF INITS AFTER INPUT BLOCK */ term = 1; while (term <= max_terms) { /* do number 2*/ array_y_init[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_norms[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_fact_1[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_pole[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_1st_rel_error[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_last_rel_error[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_type_pole[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_y[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_x[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_tmp0[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_tmp1[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_m1[term] = 0.0; term = term + 1; }/* end do number 2*/ ; ord = 1; while (ord <=2) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y_higher[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=2) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y_higher_work[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=2) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y_higher_work2[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=2) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y_set_initial[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=1) { /* do number 2*/ term = 1; while (term <= 3) { /* do number 3*/ array_poles[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=1) { /* do number 2*/ term = 1; while (term <= 3) { /* do number 3*/ array_real_pole[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=1) { /* do number 2*/ term = 1; while (term <= 3) { /* do number 3*/ array_complex_pole[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=max_terms) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_fact_2[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; /* BEGIN ARRAYS DEFINED AND INITIALIZATED */ term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_y[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_x[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_tmp0[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_tmp1[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_m1[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_const_1[term] = 0.0; term = term + 1; }/* end do number 2*/ ; array_const_1[1] = 1; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_const_0D0[term] = 0.0; term = term + 1; }/* end do number 2*/ ; array_const_0D0[1] = 0.0; term = 1; while (term <= max_terms) { /* do number 2*/ array_m1[term] = 0.0; term = term + 1; }/* end do number 2*/ ; array_m1[1] = -1.0; /* END ARRAYS DEFINED AND INITIALIZATED */ /* Initing Factorial Tables */ iiif = 0; while (iiif <= glob_max_terms) { /* do number 2*/ jjjf = 0; while (jjjf <= glob_max_terms) { /* do number 3*/ array_fact_1[iiif] = 0; array_fact_2[iiif][jjjf] = 0; jjjf = jjjf + 1; }/* end do number 3*/ ; iiif = iiif + 1; }/* end do number 2*/ ; /* Done Initing Factorial Tables */ /* TOP SECOND INPUT BLOCK */ /* BEGIN SECOND INPUT BLOCK */ /* END FIRST INPUT BLOCK */ /* BEGIN SECOND INPUT BLOCK */ x_start=-5.0; x_end=5.0; array_y_init[0 + 1] = exact_soln_y(x_start); glob_look_poles=false; glob_max_iter=1000000; glob_display_interval=0.1; glob_max_minutes=10; /* END SECOND INPUT BLOCK */ /* BEGIN OVERRIDE BLOCK */ glob_desired_digits_correct=10; glob_display_interval=0.001; glob_look_poles=true; glob_max_iter=10000000; glob_max_minutes=3; glob_subiter_method=3; /* END OVERRIDE BLOCK */ /* END SECOND INPUT BLOCK */ /* BEGIN INITS AFTER SECOND INPUT BLOCK */ glob_last_good_h = glob_h; glob_max_terms = max_terms; glob_max_sec = convfloat(60.0) * convfloat(glob_max_minutes) + convfloat(3600.0) * convfloat(glob_max_hours); if (glob_h > 0.0) { /* if number 1*/ glob_neg_h = false; glob_display_interval = omniabs(glob_display_interval); } else { glob_neg_h = true; glob_display_interval = -omniabs(glob_display_interval); }/* end if 1*/ ; chk_data(); /* AFTER INITS AFTER SECOND INPUT BLOCK */ 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; /* BEGIN OPTIMIZE CODE */ omniout_str(ALWAYS,"START of Optimize"); /* Start Series -- INITIALIZE FOR OPTIMIZE */ glob_check_sign = check_sign(x_start,x_end); glob_h = check_sign(x_start,x_end); if (glob_display_interval < glob_h) { /* if number 2*/ glob_h = glob_display_interval; }/* end if 2*/ ; if (glob_max_h < glob_h) { /* if number 2*/ glob_h = glob_max_h; }/* end if 2*/ ; found_h = -1.0; best_h = 0.0; min_value = glob_large_float; est_answer = est_size_answer(); opt_iter = 1; while ((opt_iter <= 20) && (found_h < 0.0)) { /* do number 2*/ 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 3*/ 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 3*/ ; rows = order_diff; r_order = 1; while (r_order <= rows) { /* do number 3*/ term_no = 1; while (term_no <= (rows - r_order + 1)) { /* do number 4*/ it = term_no + r_order - 1; array_y_higher[r_order][term_no] = array_y_init[it]* expt(glob_h , (term_no - 1)) / ((factorial_1(term_no - 1))); term_no = term_no + 1; }/* end do number 4*/ ; r_order = r_order + 1; }/* end do number 3*/ ; atomall(); 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,""); value3 = test_suggested_h(); omniout_float(ALWAYS,"value3",32,value3,32,""); if ((value3 < est_needed_step_err) && (found_h < 0.0)) { /* if number 2*/ best_h = glob_h; found_h = 1.0; }/* end if 2*/ ; omniout_float(ALWAYS,"best_h",32,best_h,32,""); opt_iter = opt_iter + 1; glob_h = glob_h * 0.5; }/* end do number 2*/ ; if (found_h > 0.0) { /* if number 2*/ glob_h = best_h ; } else { omniout_str(ALWAYS,"No increment to obtain desired accuracy found"); }/* end if 2*/ ; /* END OPTIMIZE CODE */ if (glob_html_log) { /* if number 2*/ html_log_file = fopen("html/entry.html","w"); }/* end if 2*/ ; /* BEGIN SOLUTION CODE */ if (found_h > 0.0) { /* if number 2*/ 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; 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; array_y_higher[r_order][term_no] = array_y_init[it]* expt(glob_h , (term_no - 1)) / ((factorial_1(term_no - 1))); term_no = term_no + 1; }/* end do number 3*/ ; r_order = r_order + 1; }/* end do number 2*/ ; 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 )) && ((convfloat(glob_clock_sec) - convfloat(glob_orig_start_sec)) < convfloat(glob_max_sec))) { /* do number 2*/ /* left paren 0001C */ if (reached_interval()) { /* if number 3*/ omniout_str(INFO," "); omniout_str(INFO,"TOP MAIN SOLVE Loop"); }/* end if 3*/ ; glob_iter = glob_iter + 1; glob_clock_sec = elapsed_time_seconds(); glob_current_iter = glob_current_iter + 1; atomall(); display_alot(current_iter); if (glob_look_poles) { /* if number 3*/ /* left paren 0004C */ check_for_pole(); }/* end if 3*/ ;/* was right paren 0004C */ if (reached_interval()) { /* if number 3*/ glob_next_display = glob_next_display + glob_display_interval; }/* end if 3*/ ; 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 = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ 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 3*/ ; /* AFTER ADJUST SUBSERIES EQ =1 */ /* BEFORE SUM SUBSERIES EQ =1 */ temp_sum = 0.0; ord = 2; calc_term = 1; /* sum_subseriesarray_y */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; 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 = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ 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 3*/ ; /* AFTER ADJUST SUBSERIES EQ =1 */ /* BEFORE SUM SUBSERIES EQ =1 */ temp_sum = 0.0; ord = 1; calc_term = 2; /* sum_subseriesarray_y */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; 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 = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ 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 3*/ ; /* AFTER ADJUST SUBSERIES EQ =1 */ /* BEFORE SUM SUBSERIES EQ =1 */ temp_sum = 0.0; ord = 1; calc_term = 1; /* sum_subseriesarray_y */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; 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 = glob_max_terms; while (term_no >= 1) { /* do number 3*/ array_y[term_no] = array_y_higher_work2[1][term_no]; ord = 1; while (ord <= order_diff) { /* do number 4*/ array_y_higher[ord][term_no] = array_y_higher_work2[ord][term_no]; ord = ord + 1; }/* end do number 4*/ ; term_no = term_no - 1; }/* end do number 3*/ ; /* END PART 2 HEVE MOVED TERMS to REGULAR Array */ ; }/* end do number 2*/ ;/* right paren 0001C */ omniout_str(ALWAYS,"Finished!"); if (glob_iter >= glob_max_iter) { /* if number 3*/ omniout_str(ALWAYS,"Maximum Iterations Reached before Solution Completed!"); }/* end if 3*/ ; if (elapsed_time_seconds() - convfloat(glob_orig_start_sec) >= convfloat(glob_max_sec )) { /* if number 3*/ omniout_str(ALWAYS,"Maximum Time Reached before Solution Completed!"); }/* end if 3*/ ; glob_clock_sec = elapsed_time_seconds(); omniout_str(INFO,"diff ( y , x , 1 ) = y;"); omniout_int(INFO,"Iterations ",32,glob_iter,4," ") ; prog_report(x_start,x_end); if (glob_html_log) { /* if number 3*/ logstart(html_log_file); logitem_str(html_log_file,"2013-01-28T12:43:01-06:00") ; logitem_str(html_log_file,"c") ; logitem_str(html_log_file,"diff0") ; logitem_str(html_log_file,"diff ( y , x , 1 ) = y;") ; 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_str(html_log_file,"16") ; logitem_good_digits(html_log_file,array_last_rel_error[1]) ; logitem_integer(html_log_file,glob_max_terms) ; logitem_float(html_log_file,array_1st_rel_error[1]) ; logitem_float(html_log_file,array_last_rel_error[1]) ; logitem_integer(html_log_file,glob_iter) ; logitem_pole(html_log_file,array_type_pole[1]) ; if (array_type_pole[1] == 1 || array_type_pole[1] == 2) { /* if number 4*/ logitem_float(html_log_file,array_pole[1]) ; logitem_float(html_log_file,array_pole[2]) ; 0; } else { logitem_str(html_log_file,"NA") ; logitem_str(html_log_file,"NA") ; 0; }/* end if 4*/ ; logitem_time(html_log_file,convfloat(glob_clock_sec)) ; if (glob_percent_done < 100.0) { /* if number 4*/ logitem_time(html_log_file,convfloat(glob_total_exp_sec)) ; 0; } else { logitem_str(html_log_file,"Done") ; 0; }/* end if 4*/ ; log_revs(html_log_file," 165 ") ; logitem_str(html_log_file,"diff0 diffeq.c") ; logitem_str(html_log_file,"diff0 c results") ; logitem_str(html_log_file,"All Tests - All Languages") ; logend(html_log_file) ; ; }/* end if 3*/ ; if (glob_html_log) { /* if number 3*/ fclose(html_log_file); }/* end if 3*/ ; ;; }/* end if 2*/ /* END OUTFILEMAIN */ } /* End Function number1 */