2 * CDE - Common Desktop Environment
4 * Copyright (c) 1993-2012, The Open Group. All rights reserved.
6 * These libraries and programs are free software; you can
7 * redistribute them and/or modify them under the terms of the GNU
8 * Lesser General Public License as published by the Free Software
9 * Foundation; either version 2 of the License, or (at your option)
12 * These libraries and programs are distributed in the hope that
13 * they will be useful, but WITHOUT ANY WARRANTY; without even the
14 * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15 * PURPOSE. See the GNU Lesser General Public License for more
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with these libraries and programs; if not, write
20 * to the Free Software Foundation, Inc., 51 Franklin Street, Fifth
21 * Floor, Boston, MA 02110-1301 USA
23 /* $XConsortium: functions.c /main/6 1996/09/25 09:36:28 mustafa $ */
26 * Contains the many of the functions (i.e. do_*) which actually do *
27 * (at least start) the calculations. *
29 * (c) Copyright 1993, 1994 Hewlett-Packard Company *
30 * (c) Copyright 1993, 1994 International Business Machines Corp. *
31 * (c) Copyright 1993, 1994 Sun Microsystems, Inc. *
32 * (c) Copyright 1993, 1994 Novell, Inc. *
42 extern char *base_str[] ; /* Strings for each base value. */
43 extern char *dtype_str[] ; /* Strings for each display mode value. */
44 extern char *mode_str[] ; /* Strings for each mode value. */
45 extern char *ttype_str[] ; /* Strings for each trig type value. */
46 extern char *vstrs[] ; /* Various strings. */
48 extern struct button buttons[] ; /* Calculator button values. */
49 extern struct button mode_buttons[] ; /* Special "mode" buttons. */
50 extern struct menu_entry menu_entries[] ; /* All the menu strings. */
52 extern Vars v ; /* Calctool variables and options. */
54 double mods[] = { 1.0, 1.0e-1, 1.0e-2, 1.0e-3, 1.0e-4,
55 1.0e-5, 1.0e-6, 1.0e-7, 1.0e-8, 1.0e-9,
56 1.0e-10, 1.0e-11, 1.0e-12, 1.0e-13, 1.0e-14,
57 1.0e-15, 1.0e-16, 1.0e-17, 1.0e-18, 1.0e-19 };
59 static void compute_i(double *target);
60 static int count_sign_changes(double *cf, int count);
64 do_accuracy(void) /* Set display accuracy. */
68 for (i = ACC_START; i <= ACC_END; i++)
69 if (v->current == menu_entries[i].val)
71 v->accuracy = char_val(v->current) ;
80 do_ascii(void) /* Convert ASCII value. */
87 do_base(void) /* Change the current base setting. */
89 if (v->current == BASE_BIN) v->base = BIN ;
90 else if (v->current == BASE_OCT) v->base = OCT ;
91 else if (v->current == BASE_DEC) v->base = DEC ;
92 else if (v->current == BASE_HEX) v->base = HEX ;
99 set_base(enum base_type base)
101 grey_buttons(v->base) ;
102 show_display(v->MPdisp_val) ;
103 set_option_menu((int) BASEITEM, (int)v->base);
105 if (v->rstate) make_registers(MEM) ;
106 if (v->frstate) make_registers(FIN) ;
110 do_business(void) /* Perform special business mode calculations. */
112 Boolean need_show = TRUE;
113 char *display_number = NULL;
114 int MPbv[MP_SIZE], MP1[MP_SIZE], MP2[MP_SIZE], MP3[MP_SIZE], MP4[MP_SIZE] ;
115 int i, len, val, accSav ;
118 if (IS_KEY(v->current, KEY_CTRM))
120 /* Cterm - FMEM1 = int (periodic interest rate).
121 * FMEM2 = Pv (present value).
122 * FMEM4 = Fv (future value).
124 * RESULT = log(FMEM4 / FMEM2) / log(1 + FMEM1)
126 if(v->MPfvals[1] == 0.0 || v->MPfvals[2] == 0.0 || v->MPfvals[4] == 0.0)
128 char *errorMsg, *tmp;
130 /* want to undraw the button first */
131 draw_button(19, 0, 4, 3, FALSE);
132 errorMsg = GETMESSAGE(5, 7, "Can't calculate 'Compound Term'\nwithout a non zero interest rate,\na non zero Present value, and\na non zero Future value.\n");
133 tmp = XtNewString(errorMsg);
140 result = log(v->MPfvals[4] / v->MPfvals[2]) / log(1.0 + (v->MPfvals[1] / 1200));
143 mpcdm(&result, v->MPdisp_val) ;
144 make_registers(FIN) ;
149 else if (IS_KEY(v->current, KEY_DDB))
152 /* Ddb - MEM0 = cost (amount paid for asset).
153 * MEM1 = salvage (value of asset at end of its life).
154 * MEM2 = life (useful life of the asset).
155 * MEM3 = period (time period for depreciation allowance).
158 * for (i = 0; i < MEM3; i++)
160 * VAL = ((MEM0 - bv) * 2) / MEM2
168 mpcmi(v->MPmvals[3], &len) ;
169 for (i = 0; i < len; i++)
171 mpsub(v->MPmvals[0], MPbv, MP1) ;
173 mpmuli(MP1, &val, MP2) ;
174 mpdiv(MP2, v->MPmvals[2], v->MPdisp_val) ;
176 mpadd(MP1, v->MPdisp_val, MPbv) ;
179 else if (IS_KEY(v->current, KEY_FV))
182 /* Fv - FMEM3 = pmt (periodic payment).
183 * FMEM1 = int (periodic interest rate).
184 * FMEM2 = Pv (present value).
185 * FMEM0 = n (number of periods).
188 if(v->MPfvals[0] == 0.0 || v->MPfvals[1] == 0.0 ||
189 (v->MPfvals[2] == 0.0 && v->MPfvals[3] == 0.0) || v->funstate == 0)
194 doerr(GETMESSAGE(5, 2, "ERROR: No Solution"));
198 /* set FV register */
199 mpcmd(v->MPdisp_val, &(v->MPfvals[4]));
203 if ((w = 1.0 + v->MPfvals[1] / (v->MPfvals[5] * 100.0)) == 1.0)
204 result = -(v->MPfvals[2] + v->MPfvals[0] * v->MPfvals[3]);
206 result = -(v->MPfvals[2] * pow(w, v->MPfvals[0]) +
207 v->MPfvals[3] * (pow(w, v->MPfvals[0]) - 1.0) *
208 pow(w, 0.0) / (w - 1.0));
209 mpcdm(&result, v->MPdisp_val) ;
211 if(strcmp(v->display, GETMESSAGE(3, 364, "Error")) != 0)
212 mpcmd(v->MPdisp_val, &(v->MPfvals[4]));
213 make_registers(FIN) ;
216 else if (IS_KEY(v->current, KEY_PMT))
219 /* Pmt - FMEM0 = prin (principal).
220 * FMEM1 = int (periodic interest rate).
223 * RESULT = FMEM0 * (FMEM1 / (1 - pow(FMEM1 + 1, -1 * FMEM2)))
226 if(v->MPfvals[0] == 0.0 || v->MPfvals[1] == 0.0 ||
227 (v->MPfvals[2] == 0.0 && v->MPfvals[4] == 0.0) || v->funstate == 0)
232 doerr(GETMESSAGE(5, 2, "ERROR: No Solution"));
236 /* set Payment register */
237 mpcmd(v->MPdisp_val, &(v->MPfvals[3]));
242 if ((w = 1.0 + v->MPfvals[1] / (v->MPfvals[5] * 100.0)) == 1.0)
243 result = -(v->MPfvals[4] + v->MPfvals[2]) / v->MPfvals[0];
245 result = -(v->MPfvals[2] * pow(w, v->MPfvals[0]) +
246 v->MPfvals[4]) * (w - 1.0) /
247 ((pow(w, v->MPfvals[0]) - 1.0) * pow(w, 0.0));
248 mpcdm(&result, v->MPdisp_val) ;
250 if(strcmp(v->display, GETMESSAGE(3, 364, "Error")) != 0)
251 mpcmd(v->MPdisp_val, &(v->MPfvals[3]));
252 make_registers(FIN) ;
255 else if (IS_KEY(v->current, KEY_PV))
258 /* Pv - FMEM0 = pmt (periodic payment).
259 * FMEM1 = int (periodic interest rate).
262 * RESULT = FMEM0 * (1 - pow(1 + FMEM1, -1 * FMEM2)) / FMEM1
265 if(v->MPfvals[0] == 0.0 || v->MPfvals[1] == 0.0 || v->MPfvals[3] == 0.0 || v->funstate == 0)
270 doerr(GETMESSAGE(5, 2, "ERROR: No Solution"));
274 /* set PV register */
275 mpcmd(v->MPdisp_val, &(v->MPfvals[2]));
279 if ((w = 1.0 + v->MPfvals[1] / (v->MPfvals[5] * 100.0)) == 1.0)
280 result = -(v->MPfvals[4] + v->MPfvals[0] * v->MPfvals[3]);
282 result = -(v->MPfvals[4] / pow(w, v->MPfvals[0]) +
283 v->MPfvals[3] * (pow(w, v->MPfvals[0]) - 1.0) *
284 pow(w, 0.0 - v->MPfvals[0]) / (w - 1.0));
285 mpcdm(&result, v->MPdisp_val) ;
287 if(strcmp(v->display, GETMESSAGE(3, 364, "Error")) != 0)
288 mpcmd(v->MPdisp_val, &(v->MPfvals[2]));
289 make_registers(FIN) ;
292 else if (IS_KEY(v->current, KEY_RATE))
294 /* Rate - MEM0 = fv (future value).
295 * MEM1 = pv (present value).
298 * RESULT = pow(MEM0 / MEM1, 1 / MEM2) - 1
300 if(v->MPfvals[0] == 0.0 || (v->MPfvals[2] == 0.0 && v->MPfvals[3] == 0.0)
301 || (v->MPfvals[3] == 0.0 && v->MPfvals[4] == 0.0)
307 doerr(GETMESSAGE(5, 2, "ERROR: No Solution"));
312 accSav = v->accuracy;
314 display_number = make_number(v->MPdisp_val, FALSE);
315 MPstr_to_num(display_number, DEC, v->MPdisp_val);
316 /* set RATE register */
317 mpcmd(v->MPdisp_val, &(v->MPfvals[1]));
318 v->accuracy = accSav;
323 compute_i(&(v->MPfvals[1]));
324 mpcdm(&(v->MPfvals[1]), v->MPdisp_val);
325 accSav = v->accuracy;
327 display_number = make_number(v->MPdisp_val, FALSE);
328 MPstr_to_num(display_number, DEC, v->MPdisp_val);
329 v->accuracy = accSav;
333 make_registers(FIN) ;
336 STRCPY(v->display, display_number);
337 set_item(DISPLAYITEM, v->display);
340 else if (IS_KEY(v->current, KEY_SLN))
343 /* Sln - MEM0 = cost (cost of the asset).
344 * MEM1 = salvage (salvage value of the asset).
345 * MEM2 = life (useful life of the asset).
347 * RESULT = (MEM0 - MEM1) / MEM2
350 mpsub(v->MPmvals[0], v->MPmvals[1], MP1) ;
351 mpdiv(MP1, v->MPmvals[2], v->MPdisp_val) ;
353 else if (IS_KEY(v->current, KEY_SYD))
356 /* Syd - MEM0 = cost (cost of the asset).
357 * MEM1 = salvage (salvage value of the asset).
358 * MEM2 = life (useful life of the asset).
359 * MEM3 = period (period for which depreciation is computed).
361 * RESULT = ((MEM0 - MEM1) * (MEM2 - MEM3 + 1)) /
362 * (MEM2 * (MEM2 + 1) / 2)
365 mpsub(v->MPmvals[2], v->MPmvals[3], MP2) ;
367 mpaddi(MP2, &val, MP3) ;
368 mpaddi(v->MPmvals[2], &val, MP2) ;
369 mpmul(v->MPmvals[2], MP2, MP4) ;
372 mpdiv(MP4, MP2, MP1) ;
373 mpdiv(MP3, MP1, MP2) ;
374 mpsub(v->MPmvals[0], v->MPmvals[1], MP1) ;
375 mpmul(MP1, MP2, v->MPdisp_val) ;
377 else if (IS_KEY(v->current, KEY_TERM))
380 /* Term - FMEM0 = pmt (periodic payment).
381 * FMEM1 = fv (future value).
382 * FMEM2 = int (periodic interest rate).
384 * RESULT = log(1 + (FMEM1 * FMEM2 / FMEM0)) / log(1 + FMEM2)
387 if(v->MPfvals[1] == 0.0 || (v->MPfvals[2] == 0.0 && v->MPfvals[4] == 0)
388 || v->MPfvals[3] == 0.0 || v->funstate == 0)
393 doerr(GETMESSAGE(5, 2, "ERROR: No Solution"));
397 /* set Term register */
398 mpcmd(v->MPdisp_val, &(v->MPfvals[0]));
402 if ((w = 1.0 + v->MPfvals[1] / (v->MPfvals[5] * 100.0)) == 1.0)
403 result = -(v->MPfvals[4] + v->MPfvals[2]) / v->MPfvals[3];
406 double wdb = pow(w, 0.0);
408 result = log((v->MPfvals[3] * wdb / (w - 1.0) - v->MPfvals[4]) /
409 (v->MPfvals[2] * pow(w, 0.0) + v->MPfvals[3] * wdb /
410 (w - 1.0))) / log(w);
413 if(strcmp(v->display, GETMESSAGE(3, 364, "Error")) != 0)
414 mpcdm(&result, v->MPdisp_val) ;
416 if(strcmp(v->display, GETMESSAGE(3, 364, "Error")) != 0)
417 mpcmd(v->MPdisp_val, &(v->MPfvals[0]));
418 make_registers(FIN) ;
421 else if (IS_KEY(v->current, KEY_PYR))
423 mpcmd(v->MPdisp_val, &(v->MPfvals[5]));
424 result = do_round(v->MPfvals[5], 0);
428 v->MPfvals[5] = result;
429 make_registers(FIN) ;
432 else if (IS_KEY(v->current, KEY_FCLR))
438 /* clear Term register */
439 mpcmd(MP1, &(v->MPfvals[0])) ;
441 /* clear %/YR register */
442 mpcmd(MP1, &(v->MPfvals[1])) ;
444 /* clear PV register */
445 mpcmd(MP1, &(v->MPfvals[2])) ;
447 /* clear Payment register */
448 mpcmd(MP1, &(v->MPfvals[3])) ;
450 /* clear FV register */
451 mpcmd(MP1, &(v->MPfvals[4])) ;
455 mpcmd(MP1, &(v->MPfvals[5])) ;
460 if (need_show == TRUE)
461 show_display(v->MPdisp_val) ;
468 do_calc(void) /* Perform arithmetic calculation and display result. */
473 /* the financial state is false - last key was not a fin. key */
476 if (!(v->opsptr && !v->show_paren)) { /* Don't do if processing parens. */
477 if (IS_KEY(v->current, KEY_EQ) && IS_KEY(v->old_cal_value, KEY_EQ)) {
479 mpstr(v->MPdisp_val, v->MPresult) ;
481 mpstr(v->MPlast_input, v->MPdisp_val) ;
486 if (!IS_KEY(v->current, KEY_EQ) && IS_KEY(v->old_cal_value, KEY_EQ))
489 if (IS_KEY(v->cur_op, KEY_COS) || /* Cos */
490 IS_KEY(v->cur_op, KEY_SIN) || /* Sin */
491 IS_KEY(v->cur_op, KEY_TAN) || /* Tan */
492 v->cur_op == '?') /* Undefined */
493 mpstr(v->MPdisp_val, v->MPresult) ;
495 else if (IS_KEY(v->cur_op, KEY_ADD)) /* Addition */
496 mpadd(v->MPresult, v->MPdisp_val, v->MPresult) ;
498 else if (IS_KEY(v->cur_op, KEY_SUB)) /* Subtraction. */
499 mpsub(v->MPresult, v->MPdisp_val, v->MPresult) ;
501 else if (v->cur_op == '*' ||
502 IS_KEY(v->cur_op, KEY_MUL)) /* Multiplication */
503 mpmul(v->MPresult, v->MPdisp_val, v->MPresult) ;
505 else if (IS_KEY(v->cur_op, KEY_DIV)) /* Division. */
506 mpdiv(v->MPresult, v->MPdisp_val, v->MPresult) ;
508 else if (IS_KEY(v->cur_op, KEY_PER)) /* % */
510 mpmul(v->MPresult, v->MPdisp_val, v->MPresult) ;
511 MPstr_to_num("0.01", DEC, MP1) ;
512 mpmul(v->MPresult, MP1, v->MPresult) ;
515 else if (IS_KEY(v->cur_op, KEY_YTOX)) /* y^x */
516 mppwr2(v->MPresult, v->MPdisp_val, v->MPresult) ;
518 else if (IS_KEY(v->cur_op, KEY_AND)) /* And */
520 mpcmd(v->MPresult, &dres) ;
521 mpcmd(v->MPdisp_val, &dval) ;
522 dres = setbool((BOOLEAN)(ibool(dres) & ibool(dval))) ;
523 mpcdm(&dres, v->MPresult) ;
526 else if (IS_KEY(v->cur_op, KEY_OR)) /* Or */
528 mpcmd(v->MPresult, &dres) ;
529 mpcmd(v->MPdisp_val, &dval) ;
530 dres = setbool((BOOLEAN)(ibool(dres) | ibool(dval))) ;
531 mpcdm(&dres, v->MPresult) ;
534 else if (IS_KEY(v->cur_op, KEY_XOR)) /* Xor */
536 mpcmd(v->MPresult, &dres) ;
537 mpcmd(v->MPdisp_val, &dval) ;
538 dres = setbool((BOOLEAN)(ibool(dres) ^ ibool(dval))) ;
539 mpcdm(&dres, v->MPresult) ;
542 else if (IS_KEY(v->cur_op, KEY_XNOR)) /* Xnor */
544 mpcmd(v->MPresult, &dres) ;
545 mpcmd(v->MPdisp_val, &dval) ;
546 dres = setbool((BOOLEAN)(~ibool(dres) ^ ibool(dval))) ;
547 mpcdm(&dres, v->MPresult) ;
550 else if (IS_KEY(v->cur_op, KEY_EQ)) /* do nothing. */ ; /* Equals */
552 show_display(v->MPresult) ;
554 if (!(IS_KEY(v->current, KEY_EQ) && IS_KEY(v->old_cal_value, KEY_EQ)))
555 mpstr(v->MPdisp_val, v->MPlast_input) ;
557 mpstr(v->MPresult, v->MPdisp_val) ;
559 v->cur_op = v->current ;
561 v->old_cal_value = v->current ;
562 v->new_input = v->key_exp = 0 ;
567 do_clear(void) /* Clear the calculator display and re-initialize. */
570 if (v->error) set_item(DISPLAYITEM, "") ;
578 if (v->current >= '0' && v->current <= '9')
580 mpstr(v->MPcon_vals[char_val(v->current)], v->MPdisp_val) ;
581 show_display(v->MPdisp_val) ;
587 do_delete(void) /* Remove the last numeric character typed. */
589 if (strlen(v->display))
590 v->display[strlen(v->display)-1] = '\0' ;
592 /* If we were entering a scientific number, and we have backspaced over
593 * the exponent sign, then this reverts to entering a fixed point number.
596 if (v->key_exp && !(strchr(v->display, '+')))
599 v->display[strlen(v->display)-1] = '\0' ;
600 set_item(OPITEM, "") ;
603 /* If we've backspaced over the numeric point, clear the pointed flag. */
605 if (v->pointed && !(strchr(v->display, '.'))) v->pointed = 0 ;
607 if(strcmp(v->display, "") == 0)
610 set_item(DISPLAYITEM, v->display) ;
611 MPstr_to_num(v->display, v->base, v->MPdisp_val) ;
616 do_exchange(void) /* Exchange display with memory register. */
618 int i, MPtemp[MP_SIZE] ;
620 for (i = MEM_START; i <= MEM_END; i++)
621 if (v->current == menu_entries[i].val)
623 mpstr(v->MPdisp_val, MPtemp) ;
624 mpstr(v->MPmvals[char_val(v->current)], v->MPdisp_val) ;
625 mpstr(MPtemp, v->MPmvals[char_val(v->current)]) ;
626 make_registers(MEM) ;
633 do_expno(void) /* Get exponential number. */
635 /* the financial state is false - last key was not a fin. key */
638 v->pointed = (strchr(v->display, '.') != NULL) ;
641 STRCPY(v->display, "1.0 +") ;
642 v->new_input = v->pointed = 1 ;
644 else if (!v->pointed)
646 STRNCAT(v->display, ". +", 3) ;
649 else if (!v->key_exp) STRNCAT(v->display, " +", 2) ;
652 v->exp_posn = strchr(v->display, '+') ;
653 set_item(DISPLAYITEM, v->display) ;
654 MPstr_to_num(v->display, v->base, v->MPdisp_val) ;
659 do_factorial(int *MPval, int *MPres) /* Calculate the factorial of MPval. */
662 int i, MPa[MP_SIZE], MP1[MP_SIZE], MP2[MP_SIZE] ;
664 /* NOTE: do_factorial, on each iteration of the loop, will attempt to
665 * convert the current result to a double. If v->error is set,
666 * then we've overflowed. This is to provide the same look&feel
669 * XXX: Needs to be improved. Shouldn't need to convert to a double in
670 * order to check this.
677 if (mpeq(MPval, MP1) && mpge(MPval, MP2)) /* Only positive integers. */
680 if (mpeq(MP1, MP2)) /* Special case for 0! */
687 if (!i) matherr((struct exception *) NULL) ;
691 mpmuli(MPa, &i, MPa) ;
693 if (v->error) break ;
697 else matherr((struct exception *) NULL) ;
703 do_frame(void) /* Exit dtcalc. */
709 do_function(void) /* Perform a user defined function. */
711 enum fcp_type scurwin ;
712 int fno, scolumn, srow ;
715 scolumn = v->column ;
716 scurwin = v->curwin ;
718 if (v->current >= '0' && v->current <= '9')
720 fno = char_val(v->current) ;
721 if(strcmp(v->fun_vals[fno], "") != 0)
722 process_str(v->fun_vals[fno], M_FUN) ;
724 v->curwin = scurwin ;
726 v->column = scolumn ;
734 int i, MP1[MP_SIZE], MP2[MP_SIZE] ;
736 /* the financial state is false - last key was not a fin. key */
739 if (IS_KEY(v->current, KEY_HYP)) /* Hyp */
741 v->hyperbolic = !v->hyperbolic ;
742 set_item(HYPITEM, (v->hyperbolic) ? vstrs[(int) V_HYP]
746 else if (IS_KEY(v->current, KEY_INV)) /* Inv */
748 v->inverse = !v->inverse ;
749 set_item(INVITEM, (v->inverse) ? vstrs[(int) V_INV]
753 else if (IS_KEY(v->current, KEY_32)) /* &32 */
755 mpcmd(v->MPdisp_val, &dval) ;
756 dval2 = ibool2(dval);
758 doerr(GETMESSAGE(5, 6, "ERR:Num too large for operation"));
761 dval = setbool((BOOLEAN)dval2) ;
762 mpcdm(&dval, v->MPdisp_val) ;
766 else if (IS_KEY(v->current, KEY_16)) /* &16 */
768 mpcmd(v->MPdisp_val, &dval) ;
769 dval2 = ibool2(dval);
771 doerr(GETMESSAGE(5, 6, "ERR:Num too large for operation"));
774 dval = setbool((BOOLEAN)(ibool(dval2) & 0xffff)) ;
775 mpcdm(&dval, v->MPdisp_val) ;
779 else if (IS_KEY(v->current, KEY_ETOX)) /* e^x */
781 mpstr(v->MPdisp_val, MP1) ;
782 mpexp(MP1, v->MPdisp_val) ;
785 else if (IS_KEY(v->current, KEY_TTOX)) /* 10^x */
789 mppwr2(MP1, v->MPdisp_val, v->MPdisp_val) ;
792 else if (IS_KEY(v->current, KEY_LN)) /* Ln */
794 mpstr(v->MPdisp_val, MP1) ;
795 mpln(MP1, v->MPdisp_val) ;
798 else if (IS_KEY(v->current, KEY_LOG)) /* Log */
800 mplog10(v->MPdisp_val, v->MPdisp_val) ;
803 else if (IS_KEY(v->current, KEY_RAND)) /* Rand */
806 mpcdm(&dval, v->MPdisp_val) ;
809 else if (IS_KEY(v->current, KEY_SQRT)) /* Sqrt */
811 mpstr(v->MPdisp_val, MP1) ;
812 mpsqrt(MP1, v->MPdisp_val) ;
815 else if (IS_KEY(v->current, KEY_NOT)) /* Not */
817 mpcmd(v->MPdisp_val, &dval) ;
818 dval = setbool((BOOLEAN)~ibool(dval)) ;
819 mpcdm(&dval, v->MPdisp_val) ;
822 else if (IS_KEY(v->current, KEY_REC)) /* 1/x */
826 mpstr(v->MPdisp_val, MP2) ;
827 mpdiv(MP1, MP2, v->MPdisp_val) ;
829 else if (IS_KEY(v->current, KEY_FACT)) /* x! */
831 do_factorial(v->MPdisp_val, MP1) ;
832 mpstr(MP1, v->MPdisp_val) ;
834 else if (IS_KEY(v->current, KEY_SQR)) /* x^2 */
836 mpstr(v->MPdisp_val, MP1) ;
837 mpmul(MP1, MP1, v->MPdisp_val) ;
840 else if (IS_KEY(v->current, KEY_CHS)) /* +/- */
844 if (*v->exp_posn == '+') *v->exp_posn = '-' ;
845 else *v->exp_posn = '+' ;
846 set_item(DISPLAYITEM, v->display) ;
847 MPstr_to_num(v->display, v->base, v->MPdisp_val) ;
852 mpneg(v->MPdisp_val, v->MPdisp_val) ;
853 mpstr(v->MPdisp_val, v->MPlast_input) ;
856 show_display(v->MPdisp_val) ;
861 do_keys(void) /* Display/undisplay the dtcalc key values. */
863 v->tstate = !v->tstate ;
868 do_mode(void) /* Set special calculator mode. */
870 if (v->current == MODE_FIN) v->modetype = FINANCIAL ;
871 else if (v->current == MODE_LOG) v->modetype = LOGICAL ;
872 else if (v->current == MODE_SCI) v->modetype = SCIENTIFIC ;
875 v->curwin = FCP_KEY ;
880 do_none(void) /* Null routine for empty buttons. */
890 static int maxvals[4] = { 1, 7, 9, 15 } ;
892 /* the financial state is false - last key was not a fin. key */
895 nextchar = v->current ;
896 n = v->current - '0' ;
897 if (v->base == HEX && v->current >= 'a' && v->current <= 'f')
899 nextchar -= 32 ; /* Convert to uppercase hex digit. */
900 n = v->current - 'a' + 10 ;
902 if (n > maxvals[(int) v->base])
910 SPRINTF(v->display, "%c", nextchar) ;
915 len = strlen(v->display) ;
916 if (len < MAX_DIGITS)
918 v->display[len] = nextchar ;
919 v->display[len+1] = '\0' ;
924 set_item(DISPLAYITEM, v->display) ;
925 MPstr_to_num(v->display, v->base, v->MPdisp_val) ;
931 do_numtype(void) /* Set number type (engineering, fixed or scientific). */
933 if (v->current == DISP_ENG) v->dtype = ENG ;
934 else if (v->current == DISP_FIX) v->dtype = FIX ;
935 else if (v->current == DISP_SCI) v->dtype = SCI ;
938 set_numtype(v->dtype);
942 set_numtype(enum num_type dtype)
945 show_display(v->MPdisp_val) ;
946 set_option_menu((int) NUMITEM, (int)v->dtype);
947 if (v->rstate) make_registers(MEM) ;
948 if (v->frstate) make_registers(FIN) ;
957 /* the financial state is false - last key was not a fin. key */
960 /* Check to see if this is the first outstanding parenthesis. If so, and
961 * their is a current operation already defined, then add the current
962 * operation to the parenthesis expression being displayed.
963 * Increment parentheses count, and add the open paren to the expression.
966 if (IS_KEY(v->current, KEY_LPAR))
968 if (v->noparens == 0)
970 /* if not in default state, put the operand between the display
971 value and the paren, else just put the paren */
974 /* there is no paren, and there is no current operand ... Let's
975 make the current operand into a "x" */
982 /* if the current op is an '=' and the result in the display is
983 zero, we want to ignore the display */
986 mpcmd(v->MPdisp_val, &tmpdb);
990 STRCPY(v->display, "") ;
991 set_item(DISPLAYITEM, v->display) ;
998 paren_disp(v->cur_op) ;
1004 paren_disp(v->cur_op) ;
1009 STRCPY(v->display, "") ;
1010 set_item(DISPLAYITEM, v->display) ;
1015 int len = strlen(v->display);
1017 if(v->display[len - 1] >= '0' && v->display[len - 1] <= '9')
1018 paren_disp(v->cur_op) ;
1021 v->pending = v->current ;
1025 /* If we haven't had any left brackets yet, and this is a right bracket,
1026 * then just ignore it.
1027 * Decrement the bracket count. If the count is zero, then process the
1028 * parenthesis expression.
1031 else if (IS_KEY(v->current, KEY_RPAR))
1033 if (!v->noparens) return ;
1038 paren_disp(v->current) ;
1040 while (*ptr != '(') ptr++ ;
1041 while (*ptr != '\0') process_parens(*ptr++) ;
1045 paren_disp(v->current) ;
1052 /* the financial state is false - last key was not a fin. key */
1055 /* Certain pending operations which are half completed, force the numeric
1056 * keypad to be reshown (assuming they already aren't).
1058 * Con, Exch, Fun, Sto, Rcl and Acc show buttons 0 - 9.
1059 * < and > show buttons 0 - f.
1064 if (IS_KEY(v->current, KEY_CON) || /* Con. */
1065 IS_KEY(v->current, KEY_EXCH) || /* Exch. */
1066 IS_KEY(v->current, KEY_FUN) || /* Fun. */
1067 IS_KEY(v->current, KEY_STO) || /* Sto. */
1068 IS_KEY(v->current, KEY_RCL) || /* Rcl. */
1069 IS_KEY(v->current, KEY_ACC)) /* Acc. */
1071 if (IS_KEY(v->current, KEY_LSFT) ||
1072 IS_KEY(v->current, KEY_RSFT))
1076 if (IS_KEY(v->pending, KEY_BASE)) do_base() ; /* Base */
1077 else if (IS_KEY(v->pending, KEY_DISP)) do_numtype() ; /* Disp */
1078 else if (IS_KEY(v->pending, KEY_TRIG)) do_trigtype() ; /* Trig */
1079 else if (IS_KEY(v->pending, KEY_CON)) do_constant() ; /* Con */
1080 else if (IS_KEY(v->pending, KEY_EXCH)) do_exchange() ; /* Exch */
1081 else if (IS_KEY(v->pending, KEY_FUN)) do_function() ; /* Fun */
1082 else if (IS_KEY(v->pending, KEY_STO) || /* Sto */
1083 IS_KEY(v->pending, KEY_RCL)) /* Rcl */
1086 if (IS_KEY(v->pending_op, KEY_ADD) ||
1087 IS_KEY(v->pending_op, KEY_SUB) ||
1088 IS_KEY(v->pending_op, KEY_MUL) ||
1089 IS_KEY(v->pending_op, KEY_DIV)) return ;
1091 else if (IS_KEY(v->pending, KEY_LSFT) || /* < */
1092 IS_KEY(v->pending, KEY_RSFT)) do_shift() ; /* > */
1093 else if (IS_KEY(v->pending, KEY_ACC)) do_accuracy() ; /* Acc */
1094 else if (IS_KEY(v->pending, KEY_MODE)) do_mode() ; /* Mode */
1095 else if (IS_KEY(v->pending, KEY_LPAR)) /* ( */
1100 else if (!v->pending)
1102 save_pending_values(v->current) ;
1103 v->pending_op = KEY_EQ ;
1107 show_display(v->MPdisp_val) ;
1108 if (v->error) set_item(OPITEM, vstrs[(int) V_CLR]) ;
1109 else set_item(OPITEM, "") ; /* Redisplay pending op. (if any). */
1113 grey_buttons(v->base) ; /* Just show numeric keys for current base. */
1118 do_point(void) /* Handle numeric point. */
1120 /* the financial state is false - last key was not a fin. key */
1127 STRCPY(v->display, ".") ;
1130 else STRNCAT(v->display, ".", 1) ;
1135 set_item(DISPLAYITEM, v->display) ;
1136 MPstr_to_num(v->display, v->base, v->MPdisp_val) ;
1145 /* the financial state is false - last key was not a fin. key */
1148 if (IS_KEY(v->current, KEY_ABS)) /* Abs */
1150 mpstr(v->MPdisp_val, MP1) ;
1151 mpabs(MP1, v->MPdisp_val) ;
1153 else if (IS_KEY(v->current, KEY_FRAC)) /* Frac */
1155 mpstr(v->MPdisp_val, MP1) ;
1156 mpcmf(MP1, v->MPdisp_val) ;
1158 else if (IS_KEY(v->current, KEY_INT)) /* Int */
1160 mpstr(v->MPdisp_val, MP1) ;
1161 mpcmim(MP1, v->MPdisp_val) ;
1163 show_display(v->MPdisp_val) ;
1168 do_shift(void) /* Perform bitwise shift on display value. */
1170 int MPtemp[MP_SIZE], shift ;
1174 shift = char_val(v->current) ;
1175 if(strcmp(v->snum, v->display) != 0)
1177 MPstr_to_num(v->display, v->base, MPtemp) ;
1178 mpcmd(MPtemp, &dval) ;
1181 mpcmd(v->MPdisp_val, &dval) ;
1182 temp = ibool(dval) ;
1184 if (IS_KEY(v->pending, KEY_LSFT)) temp = temp << shift ;
1185 else if (IS_KEY(v->pending, KEY_RSFT)) temp = temp >> shift ;
1187 dval = setbool((BOOLEAN)temp) ;
1188 mpcdm(&dval, v->MPdisp_val) ;
1189 show_display(v->MPdisp_val) ;
1190 mpstr(v->MPdisp_val, v->MPlast_input) ;
1196 do_sto_rcl(void) /* Save/restore value to/from memory register. */
1198 int i, MPn[MP_SIZE], n ;
1200 for (i = MEM_START; i <= MEM_END; i++)
1201 if (v->current == menu_entries[i].val)
1203 if (IS_KEY(v->pending, KEY_RCL)) /* Rcl */
1205 mpstr(v->MPmvals[char_val(v->current)], v->MPdisp_val) ;
1208 else if (IS_KEY(v->pending, KEY_STO)) /* Sto */
1210 n = char_val(v->current) ;
1212 if (IS_KEY(v->pending_op, KEY_ADD)) /* + */
1214 mpstr(v->MPmvals[n], MPn) ;
1215 mpadd(MPn, v->MPdisp_val, v->MPmvals[n]) ;
1217 else if (IS_KEY(v->pending_op, KEY_SUB)) /* - */
1219 mpstr(v->MPmvals[n], MPn) ;
1220 mpsub(MPn, v->MPdisp_val, v->MPmvals[n]) ;
1222 else if (IS_KEY(v->pending_op, KEY_MUL)) /* x */
1224 mpstr(v->MPmvals[n], MPn) ;
1225 mpmul(MPn, v->MPdisp_val, v->MPmvals[n]) ;
1227 else if (IS_KEY(v->pending_op, KEY_DIV)) /* / */
1229 mpstr(v->MPmvals[n], MPn) ;
1230 mpdiv(MPn, v->MPdisp_val, v->MPmvals[n]) ;
1232 else mpstr(v->MPdisp_val, v->MPmvals[n]) ;
1235 make_registers(MEM) ;
1240 if (IS_KEY(v->current, KEY_ADD) || IS_KEY(v->current, KEY_SUB) ||
1241 IS_KEY(v->current, KEY_MUL) || IS_KEY(v->current, KEY_DIV))
1242 v->pending_op = v->current ;
1247 do_trig(void) /* Perform all trigonometric functions. */
1249 int i, MPtemp[MP_SIZE], MP1[MP_SIZE], MP2[MP_SIZE] ;
1251 int MPcos[MP_SIZE], MPsin[MP_SIZE] ;
1257 if (v->ttype == DEG)
1260 mpmul(v->MPdisp_val, MP1, MP2) ;
1263 mpdiv(MP2, MP1, MPtemp) ;
1265 else if (v->ttype == GRAD)
1268 mpmul(v->MPdisp_val, MP1, MP2) ;
1271 mpdiv(MP2, MP1, MPtemp) ;
1273 else mpstr(v->MPdisp_val, MPtemp) ;
1275 else mpstr(v->MPdisp_val, MPtemp) ;
1279 if (IS_KEY(v->current, KEY_COS)) /* Cos */
1280 mpcos(MPtemp, v->MPtresults[(int) RAD]) ;
1281 else if (IS_KEY(v->current, KEY_SIN)) /* Sin */
1282 mpsin(MPtemp, v->MPtresults[(int) RAD]) ;
1283 else if (IS_KEY(v->current, KEY_TAN)) /* Tan */
1285 mpsin(MPtemp, MPsin) ;
1286 mpcos(MPtemp, MPcos) ;
1287 mpcmd(MPcos, &cval) ;
1288 if (cval == 0.0) doerr(vstrs[(int) V_ERROR]) ;
1289 mpdiv(MPsin, MPcos, v->MPtresults[(int) RAD]) ;
1294 if (IS_KEY(v->current, KEY_COS)) /* Cosh */
1295 mpcosh(MPtemp, v->MPtresults[(int) RAD]) ;
1296 else if (IS_KEY(v->current, KEY_SIN)) /* Sinh */
1297 mpsinh(MPtemp, v->MPtresults[(int) RAD]) ;
1298 else if (IS_KEY(v->current, KEY_TAN)) /* Tanh */
1299 mptanh(MPtemp, v->MPtresults[(int) RAD]) ;
1302 mpstr(v->MPtresults[(int) RAD], v->MPtresults[(int) DEG]) ;
1303 mpstr(v->MPtresults[(int) RAD], v->MPtresults[(int) GRAD]) ;
1309 if (IS_KEY(v->current, KEY_COS)) /* Acos */
1310 mpacos(v->MPdisp_val, v->MPdisp_val) ;
1311 else if (IS_KEY(v->current, KEY_SIN)) /* Asin */
1312 mpasin(v->MPdisp_val, v->MPdisp_val) ;
1313 else if (IS_KEY(v->current, KEY_TAN)) /* Atan */
1314 mpatan(v->MPdisp_val, v->MPdisp_val) ;
1318 if (IS_KEY(v->current, KEY_COS)) /* Acosh */
1319 mpacosh(v->MPdisp_val, v->MPdisp_val) ;
1320 else if (IS_KEY(v->current, KEY_SIN)) /* Asinh */
1321 mpasinh(v->MPdisp_val, v->MPdisp_val) ;
1322 else if (IS_KEY(v->current, KEY_TAN)) /* Atanh */
1323 mpatanh(v->MPdisp_val, v->MPdisp_val) ;
1330 mpmul(v->MPdisp_val, MP1, MP2) ;
1332 mpdiv(MP2, MP1, v->MPtresults[(int) DEG]) ;
1336 mpmul(v->MPdisp_val, MP1, MP2) ;
1338 mpdiv(MP2, MP1, v->MPtresults[(int) GRAD]) ;
1342 mpstr(v->MPdisp_val, v->MPtresults[(int) DEG]) ;
1343 mpstr(v->MPdisp_val, v->MPtresults[(int) GRAD]) ;
1346 mpstr(v->MPdisp_val, v->MPtresults[(int) RAD]) ;
1349 show_display(v->MPtresults[(int) v->ttype]) ;
1350 mpstr(v->MPtresults[(int) v->ttype], v->MPdisp_val) ;
1356 do_trigtype(void) /* Change the current trigonometric type. */
1358 if (v->current == TRIG_DEG) v->ttype = DEG ;
1359 else if (v->current == TRIG_GRA) v->ttype = GRAD ;
1360 else if (v->current == TRIG_RAD) v->ttype = RAD ;
1363 if (IS_KEY(v->cur_op, KEY_COS) ||
1364 IS_KEY(v->cur_op, KEY_SIN) ||
1365 IS_KEY(v->cur_op, KEY_TAN))
1367 mpstr(v->MPtresults[(int) v->ttype], v->MPdisp_val) ;
1368 show_display(v->MPtresults[(int) v->ttype]) ;
1370 set_option_menu((int) TTYPEITEM, (int)v->ttype);
1380 if (x > 68719476736.00) return(0) ;
1381 else if (x < -68719476736.00) return(0) ;
1384 while (x < 0.0) x += 4294967296.00 ;
1385 while (x >= 4294967296.00) x -= 4294967296.00 ;
1396 if (x > 9007199254740991.00 || x < -9007199254740991.00)
1402 while (x < 0.0) x += 4294967296.00 ;
1403 while (x >= 4294967296.00) x -= 4294967296.00 ;
1410 /* The following MP routines were not in the Brent FORTRAN package. They are
1411 * derived here, in terms of the existing routines.
1414 /* MP precision arc cosine.
1416 * 1. If (x < -1.0 or x > 1.0) then report DOMAIN error and return 0.0.
1418 * 2. If (x = 0.0) then acos(x) = PI/2.
1420 * 3. If (x = 1.0) then acos(x) = 0.0
1422 * 4. If (x = -1.0) then acos(x) = PI.
1424 * 5. If (0.0 < x < 1.0) then acos(x) = atan(sqrt(1-(x**2)) / x)
1426 * 6. If (-1.0 < x < 0.0) then acos(x) = atan(sqrt(1-(x**2)) / x) + PI
1430 mpacos(int *MPx, int *MPretval)
1432 int MP0[MP_SIZE], MP1[MP_SIZE], MP2[MP_SIZE] ;
1433 int MPn1[MP_SIZE], MPpi[MP_SIZE], MPy[MP_SIZE], val ;
1443 if (mpgt(MPx, MP1) || mplt(MPx, MPn1))
1445 doerr("acos DOMAIN error") ;
1446 mpstr(MP0, MPretval) ;
1448 else if (mpeq(MPx, MP0))
1451 mpdivi(MPpi, &val, MPretval) ;
1453 else if (mpeq(MPx, MP1)) mpstr(MP0, MPretval) ;
1454 else if (mpeq(MPx, MPn1)) mpstr(MPpi, MPretval) ;
1457 mpmul(MPx, MPx, MP2) ;
1458 mpsub(MP1, MP2, MP2) ;
1460 mpdiv(MP2, MPx, MP2) ;
1462 if (mpgt(MPx, MP0)) mpstr(MPy, MPretval) ;
1463 else mpadd(MPy, MPpi, MPretval) ;
1468 /* MP precision hyperbolic arc cosine.
1470 * 1. If (x < 1.0) then report DOMAIN error and return 0.0.
1472 * 2. acosh(x) = log(x + sqrt(x**2 - 1))
1476 mpacosh(int *MPx, int *MPretval)
1478 int MP1[MP_SIZE], val ;
1484 doerr("acosh DOMAIN error") ;
1486 mpcim(&val, MPretval) ;
1490 mpmul(MPx, MPx, MP1) ;
1492 mpaddi(MP1, &val, MP1) ;
1494 mpadd(MPx, MP1, MP1) ;
1495 mpln(MP1, MPretval) ;
1500 /* MP precision hyperbolic arc sine.
1502 * 1. asinh(x) = log(x + sqrt(x**2 + 1))
1506 mpasinh(int *MPx, int *MPretval)
1508 int MP1[MP_SIZE], val ;
1510 mpmul(MPx, MPx, MP1) ;
1512 mpaddi(MP1, &val, MP1) ;
1514 mpadd(MPx, MP1, MP1) ;
1515 mpln(MP1, MPretval) ;
1519 /* MP precision hyperbolic arc tangent.
1521 * 1. If (x <= -1.0 or x >= 1.0) then report a DOMAIn error and return 0.0.
1523 * 2. atanh(x) = 0.5 * log((1 + x) / (1 - x))
1527 mpatanh(int *MPx, int *MPretval)
1529 int MP0[MP_SIZE], MP1[MP_SIZE], MP2[MP_SIZE] ;
1530 int MP3[MP_SIZE], MPn1[MP_SIZE], val ;
1539 if (mpge(MPx, MP1) || mple(MPx, MPn1))
1541 doerr("atanh DOMAIN error") ;
1542 mpstr(MP0, MPretval) ;
1546 mpadd(MP1, MPx, MP2) ;
1547 mpsub(MP1, MPx, MP3) ;
1548 mpdiv(MP2, MP3, MP3) ;
1550 MPstr_to_num("0.5", DEC, MP1) ;
1551 mpmul(MP1, MP3, MPretval) ;
1556 /* MP precision common log.
1558 * 1. log10(x) = log10(e) * log(x)
1562 mplog10(int *MPx, int *MPretval)
1564 int MP1[MP_SIZE], MP2[MP_SIZE], n ;
1570 mpdiv(MP2, MP1, MPretval) ;
1575 process_parens(char current)
1578 int last_lpar ; /* Position in stack of last left paren. */
1579 int last_num ; /* Position is numeric stack to start processing. */
1581 /* Check to see if this is the first outstanding parenthesis. If so, and
1582 * their is a current operation already defined, then push the current
1583 * result on the numeric stack, and note it on the op stack, with a -1,
1584 * which has this special significance.
1585 * Zeroise current display value (in case of invalid operands inside the
1587 * Add the current pending operation to the opstack.
1588 * Increment parentheses count.
1591 if (IS_KEY(current, KEY_LPAR))
1593 if (!v->noparens && v->cur_op != '?')
1595 push_num(v->MPresult) ;
1598 mpcim(&i, v->MPdisp_val) ;
1599 push_op(v->cur_op) ;
1601 v->noparens++ ; /* Count of left brackets outstanding. */
1602 save_pending_values(current) ;
1605 /* If we haven't had any left brackets yet, and this is a right bracket,
1606 * then just ignore it.
1607 * Decrement the bracket count.
1608 * Add a equals to the op stack, to force a calculation to be performed
1609 * for two op operands. This is ignored if the preceding element of the
1610 * op stack was an immediate operation.
1611 * Work out where the preceding left bracket is in the stack, and then
1612 * process the stack from that point until this end, pushing the result
1613 * on the numeric stack, and setting the new op stack pointer appropriately.
1614 * If there are no brackets left unmatched, then clear the pending flag,
1615 * clear the stack pointers and current operation, and show the display.
1618 else if (IS_KEY(current, KEY_RPAR))
1622 last_lpar = v->opsptr - 1 ;
1623 last_num = v->numsptr ;
1624 while (!IS_KEY(v->opstack[last_lpar], KEY_LPAR))
1626 if (v->opstack[last_lpar] == -1) last_num-- ;
1629 process_stack(last_lpar + 1, last_num, v->opsptr - last_lpar - 1) ;
1635 process_stack(0, 0, v->opsptr) ;
1637 v->pending = v->opsptr = v->numsptr = 0 ;
1639 set_item(OPITEM, "") ;
1642 set_item(DISPLAYITEM, vstrs[(int) V_ERROR]) ;
1643 set_item(OPITEM, vstrs[(int) V_CLR]) ;
1644 STRCPY(v->display, vstrs[(int) V_ERROR]) ;
1648 show_display(v->MPdisp_val) ;
1649 mpstr(v->MPdisp_val, v->MPlast_input) ;
1659 push_num(int *MPval) /* Try to push value onto the numeric stack. */
1661 if (v->numsptr < 0) return ;
1662 if (v->numsptr >= MAXSTACK)
1664 STRCPY(v->display, vstrs[(int) V_NUMSTACK]) ;
1665 set_item(DISPLAYITEM, v->display) ;
1668 set_item(OPITEM, vstrs[(int) V_CLR]) ;
1672 if (v->MPnumstack[v->numsptr] == NULL)
1673 v->MPnumstack[v->numsptr] =
1674 (int *) LINT_CAST(calloc(1, sizeof(int) * MP_SIZE)) ;
1675 mpstr(MPval, v->MPnumstack[v->numsptr++]) ;
1681 push_op(int val) /* Try to push value onto the operand stack. */
1683 if (v->opsptr < 0) return ;
1684 if (v->opsptr >= MAXSTACK)
1686 STRCPY(v->display, vstrs[(int) V_OPSTACK]) ;
1687 set_item(DISPLAYITEM, v->display) ;
1689 set_item(OPITEM, vstrs[(int) V_CLR]) ;
1691 else v->opstack[v->opsptr++] = val ;
1696 save_pending_values(int val)
1701 for (n = 0; n < TITEMS; n++)
1703 if (val == buttons[n].value)
1706 v->pending_win = v->curwin ;
1707 if (v->pending_win == FCP_MODE)
1708 v->pending_mode = v->modetype ;
1718 q = p & 0x80000000 ;
1721 if (q) val += 2147483648.0 ;
1726 do_round(double result, int ndigits)
1728 char buf2[40], buffer[100];
1731 if (isnan(result)) return result;
1732 #if defined(_AIX) || defined(__aix) || defined(__sparc)
1733 temp = finite(result);
1735 return (temp > 0) ? HUGE : -HUGE;
1737 if ((temp = isinf(result))) return (temp > 0) ? HUGE : -HUGE;
1738 #endif /* _AIX or sparc*/
1740 if (ndigits >= 0 && ndigits < MAX_DIGITS)
1742 result += 0.5 * (result > 0 ? mods[ndigits] : -mods[ndigits]);
1743 result -= fmod(result, mods[ndigits]);
1746 sprintf(buf2, "%%.%dlg", MAX_DIGITS);
1747 sprintf(buffer, buf2, result);
1748 return atof(buffer);
1752 try_compute_i(double guess, double *result, int method)
1754 double sum_pos, sum_pos_prime, sum_neg, sum_neg_prime, w = guess;
1760 double term, term_prime, f, f_prime, lsp, lsn;
1762 sum_pos = sum_pos_prime = sum_neg = sum_neg_prime = 0;
1764 if (v->MPfvals[2] != 0.0)
1769 term_prime = v->MPfvals[0];
1773 term = pow(w, v->MPfvals[0]);
1774 term_prime = (v->MPfvals[0]) * pow(w, v->MPfvals[0] - 1.0);
1776 if (v->MPfvals[2] > 0.0)
1778 sum_pos += v->MPfvals[2] * term;
1779 sum_pos_prime += v->MPfvals[2] * term_prime;
1783 sum_neg -= v->MPfvals[2] * term;
1784 sum_neg_prime -= v->MPfvals[2] * term_prime;
1787 if (v->MPfvals[3] != 0.0)
1791 term = v->MPfvals[0];
1793 term_prime = v->MPfvals[0] * (v->MPfvals[0] - 1) / 2.0 +
1794 v->MPfvals[0] * (0.0);
1798 double wn = pow(w, v->MPfvals[0]);
1799 double wdb = pow(w, 0.0);
1801 term = (wn - 1.0) * wdb / (w - 1.0);
1803 term_prime = (v->MPfvals[0] * pow(w,(0.0 + v->MPfvals[0] - .01))
1804 + (wn - 1.0) * (0.0) * pow(w, (0.0 - 1.0))) /
1805 (w - 1.0) - (wn - 1.0) * wdb /
1806 ((w - 1.0) * (w - 1.0));
1809 if (v->MPfvals[3] > 0.0)
1811 sum_pos += v->MPfvals[3] * term;
1812 sum_pos_prime += v->MPfvals[3] * term_prime;
1816 sum_neg -= v->MPfvals[3] * term;
1817 sum_neg_prime -= v->MPfvals[3] * term_prime;
1820 if (v->MPfvals[4] != 0.0)
1822 if (v->MPfvals[4] > 0.0) sum_pos += v->MPfvals[4];
1823 else sum_neg -= v->MPfvals[4];
1834 f_prime = sum_pos_prime / sum_pos - sum_neg_prime / sum_neg;
1837 f = lsp / lsn - 1.0;
1838 f_prime = (lsn * sum_pos_prime / sum_pos -
1839 lsp * sum_neg_prime / sum_neg) /
1844 new_w = w - f / f_prime;
1846 #if defined(_AIX) || defined(__aix) || defined(__sparc)
1847 if (!(!isnan(new_w) && finite(new_w)))
1850 if (!(!isnan(new_w) && !isinf(new_w)))
1852 #endif /* _AIX or sparc */
1854 if (new_w == w || (w != 0.0 && fabs((new_w - w) / w) < FIN_EPSILON))
1859 if (niter++ >= MAX_FIN_ITER)
1863 *result = do_round((new_w - 1.0) * 100.0 * v->MPfvals[5], -1);
1868 compute_i(double *target)
1871 double first_period, last_period;
1878 if (first_period < 0.0 || last_period < 0.0)
1880 doerr(GETMESSAGE(5, 5, "ERROR:Invalid odd period values"));
1884 p[0] = v->MPfvals[2] + (first_period == 0.0 ? v->MPfvals[3] : 0);
1885 p[1] = v->MPfvals[3];
1886 p[2] = v->MPfvals[4] + (last_period == 0.0 ? v->MPfvals[3] : 0);
1888 nsc = count_sign_changes(p, 3);
1892 int MP1[MP_SIZE], MP2[MP_SIZE], MP3[MP_SIZE], MP4[MP_SIZE];
1897 temp = v->MPfvals[4]/v->MPfvals[2];
1901 mpcdm(&(v->MPfvals[0]), MP4);
1902 mpdiv(MP2, MP4, MP3) ;
1903 mppwr2(MP1, MP3, MP5) ;
1905 mpaddi(MP5, &val, MP1) ;
1907 mpmuli(MP1, &val, v->MPdisp_val) ;
1908 mpcmd(v->MPdisp_val, target);
1913 doerr(GETMESSAGE(5, 3, "ERROR: Multiple Solutions"));
1916 else if (v->MPfvals[0] <= 0)
1918 doerr(GETMESSAGE(5, 4, "ERROR: Term <= 0"));
1922 success = try_compute_i((double)1.0, target, 1);
1923 success = success || try_compute_i((double)1.0e-12, target, 1);
1924 success = success || try_compute_i((double)1.0, target, 2);
1925 success = success || try_compute_i((double)1.0e-12, target, 2);
1928 doerr(GETMESSAGE(5, 1, "ERROR: Computation Failed"));
1932 count_sign_changes(double *cf, int count)
1934 int i, curr_sign = 0, result = 0;
1936 for (i = 0; i < count; i++)
1938 if (cf[i] == 0.0) continue;
1942 if (cf[i] > 0.0) continue;
1946 else if (curr_sign == -1)
1948 if (cf[i] < 0.0) continue;
1954 if (cf[i] > 0.0) curr_sign = 1;
1955 else curr_sign = -1;