/* 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) (x) #define convfloat(x) (x) int glob_max_terms; long elapsed_time_seconds(); 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 exact_soln_y2(double ); double exact_soln_y1(double ); /* Top Generate Globals Definition */ double glob_log10abserr=0.0; int glob_iter=0; double glob_orig_start_sec=0.0; double glob_disp_incr=0.1; int sec_in_minute=60; int glob_good_digits=0; double glob_optimal_start=0.0; int glob_no_eqs=0; unsigned char glob_dump_analytic=false; unsigned char glob_not_yet_start_msg=true; unsigned char glob_initial_pass=true; int glob_max_opt_iter=10; int glob_curr_iter_when_opt=0; int glob_unchanged_h_cnt=0; double glob_smallish_float=0.1e-100; double glob_hmin=0.00000000001; double glob_h=0.1; int years_in_century=100; int hours_in_day=24; double glob_max_hours=0.0; double glob_relerr=0.1e-10; double glob_hmin_init=0.001; double glob_max_trunc_err=0.1e-10; int days_in_year=365; unsigned char djd_debug2=true; unsigned char glob_dump=false; double glob_max_minutes=0.0; double glob_abserr=0.1e-10; double glob_hmax=1.0; unsigned char glob_display_flag=true; unsigned char glob_warned2=false; double glob_log10_relerr=0.1e-10; double glob_large_float=9.0e100; unsigned char glob_not_yet_finished=true; int min_in_hour=60; unsigned char djd_debug=true; int glob_subiter_method=3; double glob_normmax=0.0; int glob_start=0; double glob_max_sec=10000.0; unsigned char glob_look_poles=false; double glob_clock_start_sec=0.0; double glob_optimal_expect_sec=0.1; double glob_log10relerr=0.0; unsigned char glob_warned=false; double glob_small_float=0.1e-50; int glob_max_iter=1000; double glob_almost_1=0.9990; double glob_log10normmin=0.1; int glob_current_iter=0; double glob_optimal_clock_start_sec=0.0; double glob_last_good_h=0.1; unsigned char glob_optimal_done=false; int centuries_in_millinium=10; int MAX_UNCHANGED=10; double glob_max_rel_trunc_err=0.1e-10; double glob_log10_abserr=0.1e-10; unsigned char glob_reached_optimal_h=false; double glob_clock_sec=0.0; unsigned char glob_html_log=true; double glob_percent_done=0.0; /* Bottom Generate Globals Deninition */ int sec_in_hour = (3600); int sec_in_day = (86400); int sec_in_year = (31536000); double array_fact_1[max_terms + 1]; double array_x[max_terms + 1]; double array_y2_init[max_terms + 1]; double array_m1[max_terms + 1]; double array_y2[max_terms + 1]; double array_y1[max_terms + 1]; double array_y1_init[max_terms + 1]; double array_tmp0[max_terms + 1]; double array_tmp1[max_terms + 1]; double array_tmp2[max_terms + 1]; double array_tmp3[max_terms + 1]; double array_tmp4[max_terms + 1]; double array_norms[max_terms + 1]; double array_last_rel_error[max_terms + 1]; double array_1st_rel_error[max_terms + 1]; double array_type_pole[max_terms + 1]; double array_pole[max_terms + 1]; double array_fact_2[max_terms + 1][max_terms+ 1]; double array_y2_higher[6 + 1][max_terms+ 1]; double array_poles[2 + 1][3+ 1]; double array_y2_set_initial[3 + 1][max_terms+ 1]; double array_y2_higher_work[6 + 1][max_terms+ 1]; double array_real_pole[2 + 1][3+ 1]; double array_y1_higher[2 + 1][max_terms+ 1]; double array_complex_pole[2 + 1][3+ 1]; double array_y1_higher_work[2 + 1][max_terms+ 1]; double array_y2_higher_work2[6 + 1][max_terms+ 1]; double array_y1_set_initial[3 + 1][max_terms+ 1]; double array_y1_higher_work2[2 + 1][max_terms+ 1]; double array_const_0D0[max_terms+1]; double array_const_1D0[max_terms+1]; double array_const_1[max_terms+1]; double array_const_5[max_terms+1]; void display_alot(int iter) { double abserr, analytic_val_y, ind_var, numeric_val, relerr; int term_no; /* TOP DISPLAY ALOT */ if (iter >= 0) { /* if number 1*/ ind_var = array_x[1]; omniout_float(ALWAYS,"x[1] ",33,ind_var,20," "); analytic_val_y = exact_soln_y2(ind_var); omniout_float(ALWAYS,"y2[1] (analytic) ",33,analytic_val_y,20," "); term_no = 1; numeric_val = array_y2[term_no]; abserr = omniabs(numeric_val - analytic_val_y); omniout_float(ALWAYS,"y2[1] (numeric) ",33,numeric_val,20," "); if (omniabs(analytic_val_y) != 0.0) { /* if number 2*/ relerr = abserr*100.0/omniabs(analytic_val_y); if (relerr != 0.0) { /* if number 3*/ glob_good_digits = -trunc(log10(relerr/100.0)); } else { glob_good_digits = 16; }/* end if 3*/ ; } else { relerr = -1.0 ; glob_good_digits = -1; }/* end if 2*/ ; if (glob_iter == 1) { /* if number 2*/ array_1st_rel_error[1] = relerr; } else { array_last_rel_error[1] = relerr; }/* end if 2*/ ; 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," "); ; analytic_val_y = exact_soln_y1(ind_var); omniout_float(ALWAYS,"y1[1] (analytic) ",33,analytic_val_y,20," "); term_no = 1; numeric_val = array_y1[term_no]; abserr = omniabs(numeric_val - analytic_val_y); omniout_float(ALWAYS,"y1[1] (numeric) ",33,numeric_val,20," "); if (omniabs(analytic_val_y) != 0.0) { /* if number 2*/ relerr = abserr*100.0/omniabs(analytic_val_y); if (relerr != 0.0) { /* if number 3*/ glob_good_digits = -trunc(log10(relerr/100.0)); } else { glob_good_digits = 16; }/* end if 3*/ ; } else { relerr = -1.0 ; glob_good_digits = -1; }/* end if 2*/ ; if (glob_iter == 1) { /* if number 2*/ array_1st_rel_error[2] = relerr; } else { array_last_rel_error[2] = relerr; }/* end if 2*/ ; 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," "); /* BOTTOM DISPLAY ALOT */ }/* end if 1*/ ; } 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_y2_higher[1][1]) > glob_small_float) { /* if number 1*/ tmp = omniabs(array_y2_higher[1][1]); if (tmp < glob_normmax) { /* if number 2*/ glob_normmax = tmp; }/* end if 2*/ }/* end if 1*/ ; if (omniabs(array_y1_higher[1][1]) > glob_small_float) { /* if number 1*/ tmp = omniabs(array_y1_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 */ } 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)); 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)); }/* 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 */ } void check_for_pole() { int cnt, m, n, found; double dr1, dr2, ds1, ds2, hdrc, nr1, nr2, ord_no, rad_c, rcs, rm0, rm1, rm2, rm3, rm4; /* 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 - 5 - 1; while ((m >= 10) && ((omniabs(array_y2_higher[1][m]) < glob_small_float) || (omniabs(array_y2_higher[1][m-1]) < glob_small_float) || (omniabs(array_y2_higher[1][m-2]) < glob_small_float ))) { /* do number 2*/ m = m - 1; }/* end do number 2*/ ; if (m > 10) { /* if number 1*/ rm0 = array_y2_higher[1][m]/array_y2_higher[1][m-1]; rm1 = array_y2_higher[1][m-1]/array_y2_higher[1][m-2]; hdrc = convfloat(m-1)*rm0-convfloat(m-2)*rm1; if (omniabs(hdrc) > glob_small_float) { /* if number 2*/ rcs = glob_h/hdrc; ord_no = convfloat(m-1)*rm0/hdrc - convfloat(m) + 2.0; array_real_pole[1][1] = rcs; array_real_pole[1][2] = ord_no; } else { array_real_pole[1][1] = glob_large_float; array_real_pole[1][2] = glob_large_float; }/* 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 */ /* IN RADII REAL EQ = 2 */ /* Computes radius of convergence and r_order of pole from 3 adjacent Taylor series terms. EQUATUON NUMBER 2 */ /* 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_y1_higher[1][m]) < glob_small_float) || (omniabs(array_y1_higher[1][m-1]) < glob_small_float) || (omniabs(array_y1_higher[1][m-2]) < glob_small_float ))) { /* do number 2*/ m = m - 1; }/* end do number 2*/ ; if (m > 10) { /* if number 1*/ rm0 = array_y1_higher[1][m]/array_y1_higher[1][m-1]; rm1 = array_y1_higher[1][m-1]/array_y1_higher[1][m-2]; hdrc = convfloat(m-1)*rm0-convfloat(m-2)*rm1; if (omniabs(hdrc) > glob_small_float) { /* if number 2*/ rcs = glob_h/hdrc; ord_no = convfloat(m-1)*rm0/hdrc - convfloat(m) + 2.0; array_real_pole[2][1] = rcs; array_real_pole[2][2] = ord_no; } else { array_real_pole[2][1] = glob_large_float; array_real_pole[2][2] = glob_large_float; }/* end if 2*/ } else { array_real_pole[2][1] = glob_large_float; array_real_pole[2][2] = glob_large_float; }/* end if 1*/ ; /* BOTTOM RADII REAL EQ = 2 */ /* 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 - 5 - 1; cnt = 0; while ((cnt < 5) && (n >= 10)) { /* do number 2*/ if (omniabs(array_y2_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*/ array_complex_pole[1][1] = glob_large_float; array_complex_pole[1][2] = glob_large_float; } else if ((omniabs(array_y2_higher[1][m]) >= (glob_large_float)) || (omniabs(array_y2_higher[1][m-1]) >=(glob_large_float)) || (omniabs(array_y2_higher[1][m-2]) >= (glob_large_float)) || (omniabs(array_y2_higher[1][m-3]) >= (glob_large_float)) || (omniabs(array_y2_higher[1][m-4]) >= (glob_large_float)) || (omniabs(array_y2_higher[1][m-5]) >= (glob_large_float))) { /* if number 2*/ array_complex_pole[1][1] = glob_large_float; array_complex_pole[1][2] = glob_large_float; } else { rm0 = (array_y2_higher[1][m])/(array_y2_higher[1][m-1]); rm1 = (array_y2_higher[1][m-1])/(array_y2_higher[1][m-2]); rm2 = (array_y2_higher[1][m-2])/(array_y2_higher[1][m-3]); rm3 = (array_y2_higher[1][m-3])/(array_y2_higher[1][m-4]); rm4 = (array_y2_higher[1][m-4])/(array_y2_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*/ array_complex_pole[1][1] = glob_large_float; array_complex_pole[1][2] = 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) * 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 */ /* TOP RADII COMPLEX EQ = 2 */ /* 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_y1_higher[1][n]) > glob_small_float) { /* if number 2*/ cnt = cnt + 1; } else { cnt = 0; }/* end if 2*/ ; n = n - 1; }/* end do number 2*/ ; m = n + cnt; if (m <= 10) { /* if number 2*/ array_complex_pole[2][1] = glob_large_float; array_complex_pole[2][2] = glob_large_float; } else if ((omniabs(array_y1_higher[1][m]) >= (glob_large_float)) || (omniabs(array_y1_higher[1][m-1]) >=(glob_large_float)) || (omniabs(array_y1_higher[1][m-2]) >= (glob_large_float)) || (omniabs(array_y1_higher[1][m-3]) >= (glob_large_float)) || (omniabs(array_y1_higher[1][m-4]) >= (glob_large_float)) || (omniabs(array_y1_higher[1][m-5]) >= (glob_large_float))) { /* if number 3*/ array_complex_pole[2][1] = glob_large_float; array_complex_pole[2][2] = glob_large_float; } else { rm0 = (array_y1_higher[1][m])/(array_y1_higher[1][m-1]); rm1 = (array_y1_higher[1][m-1])/(array_y1_higher[1][m-2]); rm2 = (array_y1_higher[1][m-2])/(array_y1_higher[1][m-3]); rm3 = (array_y1_higher[1][m-3])/(array_y1_higher[1][m-4]); rm4 = (array_y1_higher[1][m-4])/(array_y1_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 4*/ array_complex_pole[2][1] = glob_large_float; array_complex_pole[2][2] = glob_large_float; } else { if (omniabs(nr1*dr2 - nr2 * dr1) > glob_small_float) { /* if number 5*/ 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 6*/ if (rcs > 0.0) { /* if number 7*/ rad_c = sqrt(rcs) * glob_h; } else { rad_c = glob_large_float; }/* end if 7*/ } else { rad_c = glob_large_float; ord_no = glob_large_float; }/* end if 6*/ } else { rad_c = glob_large_float; ord_no = glob_large_float; }/* end if 5*/ }/* end if 4*/ ; array_complex_pole[2][1] = rad_c; array_complex_pole[2][2] = ord_no; }/* end if 3*/ ; /* BOTTOM RADII COMPLEX EQ = 2 */ found = false; /* TOP WHICH RADII EQ = 1 */ if ( ! found && ((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 3*/ array_poles[1][1] = array_complex_pole[1][1]; array_poles[1][2] = array_complex_pole[1][2]; found = true; array_type_pole[1] = 2; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"Complex estimate of poles used"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found && ((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] > 0.0) && ((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 3*/ array_poles[1][1] = array_real_pole[1][1]; array_poles[1][2] = array_real_pole[1][2]; found = true; array_type_pole[1] = 1; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"Real estimate of pole used"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found && (((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 3*/ array_poles[1][1] = glob_large_float; array_poles[1][2] = glob_large_float; found = true; array_type_pole[1] = 3; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"NO POLE"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found && ((array_real_pole[1][1] < array_complex_pole[1][1]) && (array_real_pole[1][1] > 0.0) && (array_real_pole[1][2] > 0.0))) { /* if number 3*/ array_poles[1][1] = array_real_pole[1][1]; array_poles[1][2] = array_real_pole[1][2]; found = true; array_type_pole[1] = 1; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"Real estimate of pole used"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found && ((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 3*/ array_poles[1][1] = array_complex_pole[1][1]; array_poles[1][2] = array_complex_pole[1][2]; array_type_pole[1] = 2; found = true; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"Complex estimate of poles used"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found ) { /* if number 3*/ array_poles[1][1] = glob_large_float; array_poles[1][2] = glob_large_float; array_type_pole[1] = 3; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"NO POLE"); }/* end if 4*/ ; }/* end if 3*/ ; /* BOTTOM WHICH RADII EQ = 1 */ found = false; /* TOP WHICH RADII EQ = 2 */ if ( ! found && ((array_real_pole[2][1] == glob_large_float) || (array_real_pole[2][2] == glob_large_float)) && ((array_complex_pole[2][1] != glob_large_float) && (array_complex_pole[2][2] != glob_large_float)) && ((array_complex_pole[2][1] > 0.0) && (array_complex_pole[2][2] > 0.0))) { /* if number 3*/ array_poles[2][1] = array_complex_pole[2][1]; array_poles[2][2] = array_complex_pole[2][2]; found = true; array_type_pole[2] = 2; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"Complex estimate of poles used"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found && ((array_real_pole[2][1] != glob_large_float) && (array_real_pole[2][2] != glob_large_float) && (array_real_pole[2][1] > 0.0) && (array_real_pole[2][2] > 0.0) && ((array_complex_pole[2][1] == glob_large_float) || (array_complex_pole[2][2] == glob_large_float) || (array_complex_pole[2][1] <= 0.0 ) || (array_complex_pole[2][2] <= 0.0)))) { /* if number 3*/ array_poles[2][1] = array_real_pole[2][1]; array_poles[2][2] = array_real_pole[2][2]; found = true; array_type_pole[2] = 1; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"Real estimate of pole used"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found && (((array_real_pole[2][1] == glob_large_float) || (array_real_pole[2][2] == glob_large_float)) && ((array_complex_pole[2][1] == glob_large_float) || (array_complex_pole[2][2] == glob_large_float)))) { /* if number 3*/ array_poles[2][1] = glob_large_float; array_poles[2][2] = glob_large_float; found = true; array_type_pole[2] = 3; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"NO POLE"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found && ((array_real_pole[2][1] < array_complex_pole[2][1]) && (array_real_pole[2][1] > 0.0) && (array_real_pole[2][2] > 0.0))) { /* if number 3*/ array_poles[2][1] = array_real_pole[2][1]; array_poles[2][2] = array_real_pole[2][2]; found = true; array_type_pole[2] = 1; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"Real estimate of pole used"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found && ((array_complex_pole[2][1] != glob_large_float) && (array_complex_pole[2][2] != glob_large_float) && (array_complex_pole[2][1] > 0.0) && (array_complex_pole[2][2] > 0.0))) { /* if number 3*/ array_poles[2][1] = array_complex_pole[2][1]; array_poles[2][2] = array_complex_pole[2][2]; array_type_pole[2] = 2; found = true; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"Complex estimate of poles used"); }/* end if 4*/ ; }/* end if 3*/ ; if ( ! found ) { /* if number 3*/ array_poles[2][1] = glob_large_float; array_poles[2][2] = glob_large_float; array_type_pole[2] = 3; if (glob_display_flag) { /* if number 4*/ omniout_str(ALWAYS,"NO POLE"); }/* end if 4*/ ; }/* end if 3*/ ; /* BOTTOM WHICH RADII EQ = 2 */ 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 3*/ array_pole[1] = array_poles[1][1]; array_pole[2] = array_poles[1][2]; }/* end if 3*/ ; /* BOTTOM WHICH RADIUS EQ = 1 */ /* TOP WHICH RADIUS EQ = 2 */ if (array_pole[1] > array_poles[2][1]) { /* if number 3*/ array_pole[1] = array_poles[2][1]; array_pole[2] = array_poles[2][2]; }/* end if 3*/ ; /* BOTTOM WHICH RADIUS EQ = 2 */ /* BOTTOM CHECK FOR POLE */ display_pole(); } void get_norms() { int iii; if ( ! glob_initial_pass) { /* if number 3*/ 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_y2[iii]) > array_norms[iii]) { /* if number 4*/ array_norms[iii] = omniabs(array_y2[iii]); }/* end if 4*/ ; iii = iii + 1; }/* end do number 2*/ ; iii = 1; while (iii <= glob_max_terms) { /* do number 2*/ if (omniabs(array_y1[iii]) > array_norms[iii]) { /* if number 4*/ array_norms[iii] = omniabs(array_y1[iii]); }/* end if 4*/ ; iii = iii + 1; }/* end do number 2*/ /* BOTTOM GET NORMS */ ; }/* end if 3*/ ; } void atomall() { int kkk, order_d, term, adj2; 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_y1[1]; /* emit pre assign xxx $eq_no = 1 i = 1 $min_hdrs = 5 */ if ( ! array_y2_set_initial[1][6]) { /* if number 1*/ if (1 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[1] * expt(glob_h , (5)) * factorial_3(0,5); array_y2[6] = temporary; array_y2_higher[1][6] = temporary; temporary = temporary / glob_h * (2.0); array_y2_higher[2][5] = temporary ; temporary = temporary / glob_h * (3.0); array_y2_higher[3][4] = temporary ; temporary = temporary / glob_h * (4.0); array_y2_higher[4][3] = temporary ; temporary = temporary / glob_h * (5.0); array_y2_higher[5][2] = temporary ; temporary = temporary / glob_h * (6.0); array_y2_higher[6][1] = temporary ; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 2; /* emit pre mult FULL FULL $eq_no = 2 i = 1 */ array_tmp3[1] = (array_m1[1] * (array_y2[1])); /* emit pre add FULL - CONST $eq_no = 2 i = 1 */ array_tmp4[1] = array_tmp3[1] + array_const_1D0[1]; /* emit pre assign xxx $eq_no = 2 i = 1 $min_hdrs = 5 */ if ( ! array_y1_set_initial[2][2]) { /* if number 1*/ if (1 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp4[1] * expt(glob_h , (1)) * factorial_3(0,1); array_y1[2] = temporary; array_y1_higher[1][2] = temporary; temporary = temporary / glob_h * (2.0); array_y1_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_y1[2]; /* emit pre assign xxx $eq_no = 1 i = 2 $min_hdrs = 5 */ if ( ! array_y2_set_initial[1][7]) { /* if number 1*/ if (2 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[2] * expt(glob_h , (5)) * factorial_3(1,6); array_y2[7] = temporary; array_y2_higher[1][7] = temporary; temporary = temporary / glob_h * (2.0); array_y2_higher[2][6] = temporary ; temporary = temporary / glob_h * (3.0); array_y2_higher[3][5] = temporary ; temporary = temporary / glob_h * (4.0); array_y2_higher[4][4] = temporary ; temporary = temporary / glob_h * (5.0); array_y2_higher[5][3] = temporary ; temporary = temporary / glob_h * (6.0); array_y2_higher[6][2] = temporary ; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 3; /* emit pre mult FULL FULL $eq_no = 2 i = 2 */ array_tmp3[2] = ats(2,array_m1,array_y2,1); /* emit pre add FULL CONST $eq_no = 2 i = 2 */ array_tmp4[2] = array_tmp3[2]; /* emit pre assign xxx $eq_no = 2 i = 2 $min_hdrs = 5 */ if ( ! array_y1_set_initial[2][3]) { /* if number 1*/ if (2 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp4[2] * expt(glob_h , (1)) * factorial_3(1,2); array_y1[3] = temporary; array_y1_higher[1][3] = temporary; temporary = temporary / glob_h * (2.0); array_y1_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_y1[3]; /* emit pre assign xxx $eq_no = 1 i = 3 $min_hdrs = 5 */ if ( ! array_y2_set_initial[1][8]) { /* if number 1*/ if (3 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[3] * expt(glob_h , (5)) * factorial_3(2,7); array_y2[8] = temporary; array_y2_higher[1][8] = temporary; temporary = temporary / glob_h * (2.0); array_y2_higher[2][7] = temporary ; temporary = temporary / glob_h * (3.0); array_y2_higher[3][6] = temporary ; temporary = temporary / glob_h * (4.0); array_y2_higher[4][5] = temporary ; temporary = temporary / glob_h * (5.0); array_y2_higher[5][4] = temporary ; temporary = temporary / glob_h * (6.0); array_y2_higher[6][3] = temporary ; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 4; /* emit pre mult FULL FULL $eq_no = 2 i = 3 */ array_tmp3[3] = ats(3,array_m1,array_y2,1); /* emit pre add FULL CONST $eq_no = 2 i = 3 */ array_tmp4[3] = array_tmp3[3]; /* emit pre assign xxx $eq_no = 2 i = 3 $min_hdrs = 5 */ if ( ! array_y1_set_initial[2][4]) { /* if number 1*/ if (3 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp4[3] * expt(glob_h , (1)) * factorial_3(2,3); array_y1[4] = temporary; array_y1_higher[1][4] = temporary; temporary = temporary / glob_h * (2.0); array_y1_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_y1[4]; /* emit pre assign xxx $eq_no = 1 i = 4 $min_hdrs = 5 */ if ( ! array_y2_set_initial[1][9]) { /* if number 1*/ if (4 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[4] * expt(glob_h , (5)) * factorial_3(3,8); array_y2[9] = temporary; array_y2_higher[1][9] = temporary; temporary = temporary / glob_h * (2.0); array_y2_higher[2][8] = temporary ; temporary = temporary / glob_h * (3.0); array_y2_higher[3][7] = temporary ; temporary = temporary / glob_h * (4.0); array_y2_higher[4][6] = temporary ; temporary = temporary / glob_h * (5.0); array_y2_higher[5][5] = temporary ; temporary = temporary / glob_h * (6.0); array_y2_higher[6][4] = temporary ; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 5; /* emit pre mult FULL FULL $eq_no = 2 i = 4 */ array_tmp3[4] = ats(4,array_m1,array_y2,1); /* emit pre add FULL CONST $eq_no = 2 i = 4 */ array_tmp4[4] = array_tmp3[4]; /* emit pre assign xxx $eq_no = 2 i = 4 $min_hdrs = 5 */ if ( ! array_y1_set_initial[2][5]) { /* if number 1*/ if (4 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp4[4] * expt(glob_h , (1)) * factorial_3(3,4); array_y1[5] = temporary; array_y1_higher[1][5] = temporary; temporary = temporary / glob_h * (2.0); array_y1_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_y1[5]; /* emit pre assign xxx $eq_no = 1 i = 5 $min_hdrs = 5 */ if ( ! array_y2_set_initial[1][10]) { /* if number 1*/ if (5 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp1[5] * expt(glob_h , (5)) * factorial_3(4,9); array_y2[10] = temporary; array_y2_higher[1][10] = temporary; temporary = temporary / glob_h * (2.0); array_y2_higher[2][9] = temporary ; temporary = temporary / glob_h * (3.0); array_y2_higher[3][8] = temporary ; temporary = temporary / glob_h * (4.0); array_y2_higher[4][7] = temporary ; temporary = temporary / glob_h * (5.0); array_y2_higher[5][6] = temporary ; temporary = temporary / glob_h * (6.0); array_y2_higher[6][5] = temporary ; }/* end if 2*/ ; }/* end if 1*/ ; kkk = 6; /* emit pre mult FULL FULL $eq_no = 2 i = 5 */ array_tmp3[5] = ats(5,array_m1,array_y2,1); /* emit pre add FULL CONST $eq_no = 2 i = 5 */ array_tmp4[5] = array_tmp3[5]; /* emit pre assign xxx $eq_no = 2 i = 5 $min_hdrs = 5 */ if ( ! array_y1_set_initial[2][6]) { /* if number 1*/ if (5 <= glob_max_terms) { /* if number 2*/ temporary = array_tmp4[5] * expt(glob_h , (1)) * factorial_3(4,5); array_y1[6] = temporary; array_y1_higher[1][6] = temporary; temporary = temporary / glob_h * (2.0); array_y1_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_y1[kkk]; /* emit assign $eq_no = 1 */ order_d = 5; if (kkk + order_d + 1 <= glob_max_terms) { /* if number 1*/ if ( ! array_y2_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_y2[kkk + order_d] = temporary; array_y2_higher[1][kkk + order_d] = temporary; term = kkk + order_d - 1; adj2 = 2; while ((adj2 <= order_d + 1) && (term >= 1)) { /* do number 2*/ temporary = temporary / glob_h * convfp(adj2); array_y2_higher[adj2][term] = temporary; adj2 = adj2 + 1; term = term - 1; }/* end do number 2*/ }/* end if 2*/ }/* end if 1*/ ; /* emit mult FULL FULL $eq_no = 2 */ array_tmp3[kkk] = ats(kkk,array_m1,array_y2,1); /* emit FULL - NOT FULL add $eq_no = 2 */ array_tmp4[kkk] = array_tmp3[kkk]; /* emit assign $eq_no = 2 */ order_d = 1; if (kkk + order_d + 1 <= glob_max_terms) { /* if number 1*/ if ( ! array_y1_set_initial[2][kkk + order_d]) { /* if number 2*/ temporary = array_tmp4[kkk] * expt(glob_h , (order_d)) / factorial_3((kkk - 1),(kkk + order_d - 1)); array_y1[kkk + order_d] = temporary; array_y1_higher[1][kkk + order_d] = temporary; term = kkk + order_d - 1; adj2 = 2; while ((adj2 <= order_d + 1) && (term >= 1)) { /* do number 2*/ temporary = temporary / glob_h * convfp(adj2); array_y1_higher[adj2][term] = temporary; adj2 = adj2 + 1; term = term - 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 ??? */ } /* BEGIN ATS LIBRARY BLOCK */ void omniout_str(int iolevel,char *str) { if (glob_iolevel >= iolevel) { printf("%s\n",str); } } void omniout_str_noeol(int iolevel,char *str) { if (glob_iolevel >= iolevel) { printf("%s",str); } } void omniout_labstr(int iolevel,char *label,char *str) { if (glob_iolevel >= iolevel) { printf("%s = %s\n",label,str); } } 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); } } } 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); } } } void omniout_float_arr(int iolevel,char *prelabel,int elemnt,int prelen,double value,int vallen,int postlabel) { if (glob_iolevel >= iolevel) { printf("%s = [ %d ] %g %s \n", prelabel,elemnt,value, postlabel); } } void dump_series(int iolevel,char *dump_label,char *series_name, double *array_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,array_series[i]); i++; } } } void cs_info(int iolevel,char *str) { } // global glob_iolevel,glob_correct_start_flag,glob_h,glob_reached_optimal_h; // if (glob_iolevel >= iolevel) then // print("cs_info " , str , " glob_correct_start_flag = " , glob_correct_start_flag , "glob_h := " , glob_h , "glob_reached_optimal_h := " , glob_reached_optimal_h) // fi; 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 / sec_in_year); sec_temp = (secs_in % sec_in_year); days_int = (sec_temp / sec_in_day) ; sec_temp = (sec_temp % sec_in_day) ; hours_int = (sec_temp / sec_in_hour); sec_temp = (sec_temp % sec_in_hour); minutes_int = (sec_temp / sec_in_minute); sec_int = (sec_temp % 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,""); } 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 / sec_in_year); sec_temp = (secs_in % sec_in_year); days_int = (sec_temp / sec_in_day) ; sec_temp = (sec_temp % sec_in_day) ; hours_int = (sec_temp / sec_in_hour); sec_temp = (sec_temp % sec_in_hour); minutes_int = (sec_temp / sec_in_minute); sec_int = (sec_temp % 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*/ } double ats( int mmm_ats,double *array_a,double *array_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 + array_a[iii_ats]*array_b[lll_ats]; iii_ats = iii_ats + 1; }/* end do number 1*/ }/* end if 7*/ ; return ret_ats; } double att( int mmm_att,double *array_aa,double *array_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 + array_aa[iii_att]*array_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; } 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*/ } void logditto(FILE * file){ fprintf(file,""); fprintf(file,"ditto"); fprintf(file,""); } void logitem_integer(FILE * file,int n){ fprintf(file,""); fprintf(file,"%d",n); fprintf(file,""); } void logitem_str(FILE * file,char * str){ fprintf(file,""); fprintf(file,str); fprintf(file,""); } 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.0) { /* if number 8*/ good_digits = -trunc(log10(rel_error/100.0)); 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,""); } void log_revs(FILE * file,char * revs){ fprintf(file,revs); } void logitem_float(FILE * file,double x){ fprintf(file,""); fprintf(file,"%g",x); fprintf(file,""); } 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,""); } void logstart(FILE * file){ fprintf(file,""); } void logend(FILE * file){ fprintf(file,"\n"); } 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*/ } 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; } 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; } 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*/ } 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; } 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; } long elapsed_time_seconds() {; struct timeval t; struct timezone tz; gettimeofday(&t,&tz); return (t.tv_sec + t.tv_usec/1000); } double expt(double x,double y) { return pow(x,y); } double arcsin(double x) { return asin(x); } double arccos(double x) { return acos(x); } double arctan(double x) { return atan(x); } double omniabs(double x) { return fabs(x); } double ln(double x) { return log(x); } double Si(double x) { return (0.0); } double Ci(double x) { return (0.0); } /* END ATS LIBRARY BLOCK */ /* BEGIN USER DEF BLOCK */ /* BEGIN USER DEF BLOCK */ double exact_soln_y1 (double x) { return(1.0 + cos(x)); } double exact_soln_y2 (double x) { return(1.0 + sin(x)); } double exact_soln_y2p (double x) { return( cos(x)); } double exact_soln_y2pp (double x) { return( -sin(x)); } double exact_soln_y2ppp (double x) { return( -cos(x)); } double exact_soln_y2pppp (double x) { return( sin(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; double temp_sum; double x_start;double x_end; double log10norm, 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 = 2; 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/mtest7postode.ode#################"); omniout_str(ALWAYS,"diff ( y2 , x , 5 ) = y1 ;"); omniout_str(ALWAYS,"diff ( y1 , x , 1 ) = m1 * y2 + 1.0;"); 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 = 0.0;"); omniout_str(ALWAYS,"x_end = 5.0;"); omniout_str(ALWAYS,"array_y1_init[0 + 1] = exact_soln_y1(x_start);"); omniout_str(ALWAYS,"array_y2_init[0 + 1] = exact_soln_y2(x_start);"); omniout_str(ALWAYS,"array_y2_init[1 + 1] = exact_soln_y2p(x_start);"); omniout_str(ALWAYS,"array_y2_init[2 + 1] = exact_soln_y2pp(x_start);"); omniout_str(ALWAYS,"array_y2_init[3 + 1] = exact_soln_y2ppp(x_start);"); omniout_str(ALWAYS,"array_y2_init[4 + 1] = exact_soln_y2pppp(x_start);"); omniout_str(ALWAYS,"glob_h = 0.00001;"); omniout_str(ALWAYS,"glob_look_poles = true;"); omniout_str(ALWAYS,"glob_max_iter = 20;"); omniout_str(ALWAYS,"/* END SECOND INPUT BLOCK */"); omniout_str(ALWAYS,"/* BEGIN OVERRIDE BLOCK */"); omniout_str(ALWAYS,"glob_h = 0.00001 ;"); omniout_str(ALWAYS,"glob_look_poles = true;"); omniout_str(ALWAYS,"glob_max_iter = 100;"); omniout_str(ALWAYS,"glob_max_minutes = 1;"); omniout_str(ALWAYS,"/* END OVERRIDE BLOCK */"); omniout_str(ALWAYS,"!"); omniout_str(ALWAYS,"/* BEGIN USER DEF BLOCK */"); omniout_str(ALWAYS,"double exact_soln_y1 (double x) { "); omniout_str(ALWAYS,"return(1.0 + cos(x));"); omniout_str(ALWAYS,"}"); omniout_str(ALWAYS,"double exact_soln_y2 (double x) { "); omniout_str(ALWAYS,"return(1.0 + sin(x));"); omniout_str(ALWAYS,"}"); omniout_str(ALWAYS,"double exact_soln_y2p (double x) { "); omniout_str(ALWAYS,"return( cos(x));"); omniout_str(ALWAYS,"}"); omniout_str(ALWAYS,"double exact_soln_y2pp (double x) { "); omniout_str(ALWAYS,"return( -sin(x));"); omniout_str(ALWAYS,"}"); omniout_str(ALWAYS,"double exact_soln_y2ppp (double x) { "); omniout_str(ALWAYS,"return( -cos(x));"); omniout_str(ALWAYS,"}"); omniout_str(ALWAYS,"double exact_soln_y2pppp (double x) { "); omniout_str(ALWAYS,"return( sin(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; glob_log10_abserr = -8.0; glob_log10_relerr = -8.0; glob_hmax = 0.01; /* 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_fact_1[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_y2_init[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*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_y2[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_y1[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_y1_init[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_tmp2[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_tmp3[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms) { /* do number 2*/ array_tmp4[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_last_rel_error[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_type_pole[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*/ ; 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*/ ; ord = 1; while (ord <=6) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y2_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 <= 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 <=3) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y2_set_initial[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=6) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y2_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 <= 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 <=2) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y1_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 <= 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 <=2) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y1_higher_work[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=6) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y2_higher_work2[ord][term] = 0.0; term = term + 1; }/* end do number 3*/ ; ord = ord + 1; }/* end do number 2*/ ; ord = 1; while (ord <=3) { /* do number 2*/ term = 1; while (term <= max_terms) { /* do number 3*/ array_y1_set_initial[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_y1_higher_work2[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_x[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_y1[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_y2[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_tmp4[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_tmp3[term] = 0.0; term = term + 1; }/* end do number 2*/ ; term = 1; while (term <= max_terms + 1) { /* do number 2*/ array_tmp2[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_tmp0[term] = 0.0; term = term + 1; }/* end do number 2*/ ; 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 + 1) { /* do number 2*/ array_const_1D0[term] = 0.0; term = term + 1; }/* end do number 2*/ ; array_const_1D0[1] = 1.0; 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_5[term] = 0.0; term = term + 1; }/* end do number 2*/ ; array_const_5[1] = 5; 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 = 0.0; x_end = 5.0; array_y1_init[0 + 1] = exact_soln_y1(x_start); array_y2_init[0 + 1] = exact_soln_y2(x_start); array_y2_init[1 + 1] = exact_soln_y2p(x_start); array_y2_init[2 + 1] = exact_soln_y2pp(x_start); array_y2_init[3 + 1] = exact_soln_y2ppp(x_start); array_y2_init[4 + 1] = exact_soln_y2pppp(x_start); glob_h = 0.00001; glob_look_poles = true; glob_max_iter = 20; /* END SECOND INPUT BLOCK */ /* BEGIN OVERRIDE BLOCK */ glob_h = 0.00001 ; glob_look_poles = true; glob_max_iter = 100; glob_max_minutes = 1; /* 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); glob_abserr = expt(10.0 , (glob_log10_abserr)); glob_relerr = expt(10.0 , (glob_log10_relerr)); chk_data(); /* AFTER INITS AFTER SECOND INPUT BLOCK */ array_y2_set_initial[1][1] = true; array_y2_set_initial[1][2] = true; array_y2_set_initial[1][3] = true; array_y2_set_initial[1][4] = true; array_y2_set_initial[1][5] = true; array_y2_set_initial[1][6] = false; array_y2_set_initial[1][7] = false; array_y2_set_initial[1][8] = false; array_y2_set_initial[1][9] = false; array_y2_set_initial[1][10] = false; array_y2_set_initial[1][11] = false; array_y2_set_initial[1][12] = false; array_y2_set_initial[1][13] = false; array_y2_set_initial[1][14] = false; array_y2_set_initial[1][15] = false; array_y2_set_initial[1][16] = false; array_y2_set_initial[1][17] = false; array_y2_set_initial[1][18] = false; array_y2_set_initial[1][19] = false; array_y2_set_initial[1][20] = false; array_y2_set_initial[1][21] = false; array_y2_set_initial[1][22] = false; array_y2_set_initial[1][23] = false; array_y2_set_initial[1][24] = false; array_y2_set_initial[1][25] = false; array_y2_set_initial[1][26] = false; array_y2_set_initial[1][27] = false; array_y2_set_initial[1][28] = false; array_y2_set_initial[1][29] = false; array_y2_set_initial[1][30] = false; array_y1_set_initial[2][1] = true; array_y1_set_initial[2][2] = false; array_y1_set_initial[2][3] = false; array_y1_set_initial[2][4] = false; array_y1_set_initial[2][5] = false; array_y1_set_initial[2][6] = false; array_y1_set_initial[2][7] = false; array_y1_set_initial[2][8] = false; array_y1_set_initial[2][9] = false; array_y1_set_initial[2][10] = false; array_y1_set_initial[2][11] = false; array_y1_set_initial[2][12] = false; array_y1_set_initial[2][13] = false; array_y1_set_initial[2][14] = false; array_y1_set_initial[2][15] = false; array_y1_set_initial[2][16] = false; array_y1_set_initial[2][17] = false; array_y1_set_initial[2][18] = false; array_y1_set_initial[2][19] = false; array_y1_set_initial[2][20] = false; array_y1_set_initial[2][21] = false; array_y1_set_initial[2][22] = false; array_y1_set_initial[2][23] = false; array_y1_set_initial[2][24] = false; array_y1_set_initial[2][25] = false; array_y1_set_initial[2][26] = false; array_y1_set_initial[2][27] = false; array_y1_set_initial[2][28] = false; array_y1_set_initial[2][29] = false; array_y1_set_initial[2][30] = false; if (glob_html_log) { /* if number 3*/ html_log_file = fopen("html/entry.html","w"); }/* end if 3*/ ; /* BEGIN SOLUTION CODE */ omniout_str(ALWAYS,"START of Soultion"); /* Start Series -- INITIALIZE FOR SOLUTION */ array_x[1] = x_start; array_x[2] = glob_h; order_diff = 5; /* Start Series array_y2 */ term_no = 1; while (term_no <= order_diff) { /* do number 2*/ array_y2[term_no] = array_y2_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_y2_higher[r_order][term_no] = array_y2_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*/ ; order_diff = 1; /* Start Series array_y1 */ term_no = 1; while (term_no <= order_diff) { /* do number 2*/ array_y1[term_no] = array_y1_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_y1_higher[r_order][term_no] = array_y1_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(); if (omniabs(array_y2_higher[1][1]) > glob_small_float) { /* if number 3*/ tmp = omniabs(array_y2_higher[1][1]); log10norm = (log10(tmp)); if (log10norm < glob_log10normmin) { /* if number 4*/ glob_log10normmin = log10norm; }/* end if 4*/ }/* end if 3*/ ; display_alot(current_iter) ; if (omniabs(array_y1_higher[1][1]) > glob_small_float) { /* if number 3*/ tmp = omniabs(array_y1_higher[1][1]); log10norm = (log10(tmp)); if (log10norm < glob_log10normmin) { /* if number 4*/ glob_log10normmin = log10norm; }/* end if 4*/ }/* end if 3*/ ; display_alot(current_iter) ; 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) && (array_x[1] <= x_end ) && ((convfloat(glob_clock_sec) - convfloat(glob_orig_start_sec)) < convfloat(glob_max_sec))) { /* do number 2*/ /* left paren 0001C */ omniout_str(INFO," "); omniout_str(INFO,"TOP MAIN SOLVE Loop"); glob_iter = glob_iter + 1; glob_clock_sec = elapsed_time_seconds(); glob_current_iter = glob_current_iter + 1; if (glob_subiter_method == 1 ) { /* if number 3*/ atomall(); } else if (glob_subiter_method == 2 ) { /* if number 4*/ subiter = 1; while (subiter <= 6) { /* do number 3*/ atomall(); subiter = subiter + 1; }/* end do number 3*/ ; } else { subiter = 1; while (subiter <= 6 + glob_max_terms) { /* do number 3*/ atomall(); subiter = subiter + 1; }/* end do number 3*/ ; }/* end if 4*/ ; if (glob_look_poles) { /* if number 4*/ /* left paren 0004C */ check_for_pole(); }/* end if 4*/ ;/* was right paren 0004C */ array_x[1] = array_x[1] + glob_h; array_x[2] = glob_h; /* Jump Series array_y2 */ order_diff = 5; /* START PART 1 SUM AND ADJUST */ /* START SUM AND ADJUST EQ =1 */ /* sum_and_adjust array_y2 */ /* BEFORE ADJUST SUBSERIES EQ =1 */ ord = 6; calc_term = 1; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[6][iii] = array_y2_higher[6][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 = 6; calc_term = 1; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 5; calc_term = 2; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[5][iii] = array_y2_higher[5][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 = 5; calc_term = 2; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 5; calc_term = 1; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[5][iii] = array_y2_higher[5][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 = 5; calc_term = 1; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 4; calc_term = 3; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[4][iii] = array_y2_higher[4][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 = 4; calc_term = 3; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 4; calc_term = 2; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[4][iii] = array_y2_higher[4][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 = 4; calc_term = 2; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 4; calc_term = 1; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[4][iii] = array_y2_higher[4][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 = 4; calc_term = 1; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 3; calc_term = 4; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[3][iii] = array_y2_higher[3][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 = 3; calc_term = 4; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 3; calc_term = 3; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[3][iii] = array_y2_higher[3][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 = 3; calc_term = 3; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 3; calc_term = 2; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[3][iii] = array_y2_higher[3][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 = 3; calc_term = 2; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 3; calc_term = 1; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[3][iii] = array_y2_higher[3][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 = 3; calc_term = 1; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 2; calc_term = 5; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[2][iii] = array_y2_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 = 5; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 2; calc_term = 4; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[2][iii] = array_y2_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 = 4; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 2; calc_term = 3; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[2][iii] = array_y2_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 = 3; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 2; calc_term = 2; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[2][iii] = array_y2_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 = 2; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 2; calc_term = 1; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[2][iii] = array_y2_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_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 6; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[1][iii] = array_y2_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 = 6; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 5; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[1][iii] = array_y2_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 = 5; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 4; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[1][iii] = array_y2_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 = 4; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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 = 3; /* adjust_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[1][iii] = array_y2_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 = 3; /* sum_subseriesarray_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[1][iii] = array_y2_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_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y2_higher_work[1][iii] = array_y2_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_y2 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y2_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y2_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_y2[term_no] = array_y2_higher_work2[1][term_no]; ord = 1; while (ord <= order_diff) { /* do number 4*/ array_y2_higher[ord][term_no] = array_y2_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 */ /* Jump Series array_y1 */ order_diff = 1; /* START PART 1 SUM AND ADJUST */ /* START SUM AND ADJUST EQ =2 */ /* sum_and_adjust array_y1 */ /* BEFORE ADJUST SUBSERIES EQ =2 */ ord = 2; calc_term = 1; /* adjust_subseriesarray_y1 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y1_higher_work[2][iii] = array_y1_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 =2 */ /* BEFORE SUM SUBSERIES EQ =2 */ temp_sum = 0.0; ord = 2; calc_term = 1; /* sum_subseriesarray_y1 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y1_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y1_higher_work2[ord][calc_term] = temp_sum * expt(glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1)); /* AFTER SUM SUBSERIES EQ =2 */ /* BEFORE ADJUST SUBSERIES EQ =2 */ ord = 1; calc_term = 2; /* adjust_subseriesarray_y1 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y1_higher_work[1][iii] = array_y1_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 =2 */ /* BEFORE SUM SUBSERIES EQ =2 */ temp_sum = 0.0; ord = 1; calc_term = 2; /* sum_subseriesarray_y1 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y1_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y1_higher_work2[ord][calc_term] = temp_sum * expt(glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1)); /* AFTER SUM SUBSERIES EQ =2 */ /* BEFORE ADJUST SUBSERIES EQ =2 */ ord = 1; calc_term = 1; /* adjust_subseriesarray_y1 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ array_y1_higher_work[1][iii] = array_y1_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 =2 */ /* BEFORE SUM SUBSERIES EQ =2 */ temp_sum = 0.0; ord = 1; calc_term = 1; /* sum_subseriesarray_y1 */ iii = glob_max_terms; while (iii >= calc_term) { /* do number 3*/ temp_sum = temp_sum + array_y1_higher_work[ord][iii]; iii = iii - 1; }/* end do number 3*/ ; array_y1_higher_work2[ord][calc_term] = temp_sum * expt(glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1)); /* AFTER SUM SUBSERIES EQ =2 */ /* END SUM AND ADJUST EQ =2 */ /* END PART 1 */ /* START PART 2 MOVE TERMS to REGULAR Array */ term_no = glob_max_terms; while (term_no >= 1) { /* do number 3*/ array_y1[term_no] = array_y1_higher_work2[1][term_no]; ord = 1; while (ord <= order_diff) { /* do number 4*/ array_y1_higher[ord][term_no] = array_y1_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 */ display_alot(current_iter) ; }/* end do number 2*/ ;/* right paren 0001C */ omniout_str(ALWAYS,"Finished!"); if (glob_iter >= glob_max_iter) { /* if number 4*/ omniout_str(ALWAYS,"Maximum Iterations Reached before Solution Completed!"); }/* end if 4*/ ; if (elapsed_time_seconds() - convfloat(glob_orig_start_sec) >= convfloat(glob_max_sec )) { /* if number 4*/ omniout_str(ALWAYS,"Maximum Time Reached before Solution Completed!"); }/* end if 4*/ ; glob_clock_sec = elapsed_time_seconds(); omniout_str(INFO,"diff ( y2 , x , 5 ) = y1 ;"); omniout_str(INFO,"diff ( y1 , x , 1 ) = m1 * y2 + 1.0;"); omniout_int(INFO,"Iterations ",32,glob_iter,4," ") ; prog_report(x_start,x_end); if (glob_html_log) { /* if number 4*/ logstart(html_log_file); logitem_str(html_log_file,"2012-09-02T22:01:06-05:00") ; logitem_str(html_log_file,"c") ; logitem_str(html_log_file,"mtest7") ; logitem_str(html_log_file,"diff ( y2 , x , 5 ) = y1 ;") ; 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 5*/ 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 5*/ ; logitem_time(html_log_file,convfloat(glob_clock_sec)) ; if (glob_percent_done < 100.0) { /* if number 5*/ logitem_time(html_log_file,convfloat(glob_optimal_expect_sec)) ; 0; } else { logitem_str(html_log_file,"Done") ; 0; }/* end if 5*/ ; log_revs(html_log_file," 126 ") ; logitem_str(html_log_file,"mtest7 diffeq.c") ; logitem_str(html_log_file,"mtest7 c results") ; logitem_str(html_log_file,"c c++ Maple and Maxima") ; logend(html_log_file) ; logditto(html_log_file) ; logditto(html_log_file) ; logditto(html_log_file) ; logitem_str(html_log_file,"diff ( y1 , x , 1 ) = m1 * y2 + 1.0;") ; logditto(html_log_file) ; logditto(html_log_file) ; logditto(html_log_file) ; logditto(html_log_file) ; logditto(html_log_file) ; logitem_good_digits(html_log_file,array_last_rel_error[2]) ; logditto(html_log_file) ; logitem_float(html_log_file,array_1st_rel_error[2]) ; logitem_float(html_log_file,array_last_rel_error[2]) ; logditto(html_log_file) ; logitem_pole(html_log_file,array_type_pole[2]) ; if (array_type_pole[2] == 1 || array_type_pole[2] == 2) { /* if number 5*/ 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 5*/ ; logditto(html_log_file) ; if (glob_percent_done < 100.0) { /* if number 5*/ logditto(html_log_file) ; 0; } else { logditto(html_log_file) ; 0; }/* end if 5*/ ; logditto(html_log_file); ; logditto(html_log_file) ; logditto(html_log_file) ; logditto(html_log_file) ; logend(html_log_file) ; ; }/* end if 4*/ ; if (glob_html_log) { /* if number 4*/ fclose(html_log_file); }/* end if 4*/ ; ;; /* END OUTFILEMAIN */ }