1(kill(all),
2  %start : absolute_real_time(),
3  algexact : true,
4  current_fpprec : fpprec,
5  fpprec : 16,
6  approx_equal(x,y,epsilon) := block(
7    [vx:setify(listofvars(x)),vy:setify(listofvars(y)),prederror:true],
8    if vx#vy then return(false),
9    x : sort(flatten(listify(x))), y : sort(flatten(listify(y))),
10    every(map(lambda([u,v],is(abs(rectform(rhs(u)-rhs(v))) < epsilon)),x,y))),
11  load("mnewton"),
12  0);
130$
14
15approx_equal(
16  mnewton(x-1.0,x,0),
17  [[x = 1.0]],
18  newtonepsilon
19  )$
20true $
21
22approx_equal(
23  mnewton(sin(x),x,1.0),
24  [[x = 0.0]],
25  newtonepsilon
26  )$
27true $
28
29
30/* from manual */
31approx_equal(
32  mnewton([x1+3*log(x1)-x2^2, 2*x1^2-x1*x2-5*x1+1], [x1, x2], [5, 5]),
33  [[x1 = 3.756834008012769, x2 = 2.779849592817897]],
34  newtonepsilon
35  )$
36true $
37
38/* from manual */
39approx_equal(
40  mnewton([2*a^a-5],[a],[1]),
41  [[a = 1.70927556786144]],
42  newtonepsilon
43  )$
44true $
45
46/* from manual */
47approx_equal(
48  mnewton([2*3^u-v/u-5, u+2^v-4], [u, v], [2, 2]),
49  [[u = 1.066618389595407, v = 1.552564766841786]],
50  newtonepsilon
51  )$
52true $
53
54/* see https://sourceforge.net/tracker/index.php?func=detail&aid=1725549&group_id=4933&atid=104933 */
55approx_equal(
56  mnewton([2*(x1-1.1234),3*(x2-2.2345)], [x1, x2], [1,2]),
57  [[x1=1.1234,x2=2.2345]],
58  newtonepsilon
59  )$
60true $
61
62
63block([newtonepsilon:bfloat(newtonepsilon)],
64  approx_equal(
65    mnewton(cos(x)-%i,x,1b0),
66    [[x = 1.570796326794897b0 - 8.81373587019543b-1*%i]],
67    newtonepsilon
68    ))$
69true $
70
71/* see http://www.math.utexas.edu/pipermail/maxima/2010/021238.html and
72https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2994196&group_id=4933 */
73block([f:diff(acos(x^2+8.363+267)/(sqrt(x^2-4.29*x+1042)*(x^2+12.8102*x+1177)),x)],
74  approx_equal(
75    mnewton(f,x,-5.0),
76    [[x = 5.69311809542034e-9*%i - 5.626466372932345]],
77    newtonepsilon
78    ))$
79true $
80
81block([fpprec:100, newtonepsilon:bfloat(newtonepsilon)*1b-8,f:diff(acos(x^2+8.363+267)/(sqrt(x^2-4.29*x+1042)*(x^2+12.8102*x+1177)),x)],
82  approx_equal(
83    mnewton(f,x,-5.0),
84    [[x = - 5.626466372928149b0]],
85    newtonepsilon
86    ))$
87true $
88
89block([fpprec:100,newtonepsilon:bfloat(10^(-fpprec/2)),f:diff(acos(x^2+8.363+267)/(sqrt(x^2-4.29*x+1042)*(x^2+12.8102*x+1177)),x)],
90  approx_equal(
91    mnewton(f,x,-5.0),
92    [[x = - 5.62646637292814898636208310765236683228303380080106742342217423415231803837099458924175436056796413b0]],
93    newtonepsilon
94    ))$
95true $
96
97/* see https://sourceforge.net/tracker/index.php?func=detail&aid=1666011&group_id=4933&atid=104933 */
98approx_equal(
99  mnewton([x^y - %i,x+y-2],[x,y],[1,1.5]),
100  [[x=0.563519172421*%i+0.451201761049,y=1.5487982389505-0.563519172421*%i]],
101  newtonepsilon
102  )$
103true $
104
105/* see https://sourceforge.net/tracker/index.php?func=detail&aid=1662070&group_id=4933&atid=104933 */
106approx_equal(
107  mnewton([x[1] + x[2]-1, x[1] - x[2] - 10],[x[1],x[2]],[1,2]),
108  [[x[1] = 5.5,x[2] = -4.5]],
109  newtonepsilon
110  )$
111true $
112
113/* A test for a command that failed in 2010, see
114 https://def.fe.up.pt/pipermail/maxima-discuss/2010/033075.html and
115 https://def.fe.up.pt/pipermail/maxima-discuss/2010/033080.html
116
117 As the maxima-discuss thread "COERCE-FLOAT-FUN in complexfield"
118 in Mid-Jan 2019 shows the remedy for this made the following
119 command slow and space-hungry in sbcl:
120
121     w[i,j] := random (1.0) + %i * random (1.0);
122     showtime : true$
123     M : genmatrix (w, 100, 100)$
124     lu_factor (M, complexfield)$
125     lu_factor (M, generalring)$
126
127 The reason for that was using coerce-float-fun for each matrix
128 entry which created lisp functions sbcl automatically compiled.
129*/
130approx_equal(
131  mnewton(diff(acos(x^2+8.363+267)/(sqrt(x^2-4.29*x+1042)*(x^2+12.8102*x+1177)),x),x,5),
132  [[x=-1.012828847267113*10^-12*%i-5.626466372928149]],
133  newtonepsilon
134);
135true $
136
137kill(all)$
138done $
139
140(reset(algexact),0);
1410$