## 9122 Reputation

5 years, 233 days

## workaround...

Statistics is a modern package, full of modules and generated modules, not easy to debug.
For a workaround, as you noticed, it's possible to use (for the moment)

`HR:=(X, t) -> (PDF(X, t)/(1 - CDF(X, t)));`

HR(1/2*X1 + 1/2*X2, 0.4);
2.352941176

## colon, print...

Use Tabulate(DF):  instead of  print(Tabulate(DF));
(notice the colon).

## Almost correct...

It is almost correct, but why don't you read at least a tutorial first?

```Crunch := proc(x::positive)
if x >= 100 then
return x/100
else
return x + Crunch(10*x)
fi
end proc;```

(if you allow x <=0, it results an infinite recursion).

## satisfies...

You can use the (pseudo)type satisfies

```restart;
test:=table([1=[1,2,3],2=[6,5,4]]):
mytype:=satisfies(e -> e::'table', e -> [indices(e,nolist)]::list(integer),  e -> [entries(e,nolist)]::list(list(integer))):
type(test, mytype);
```

true

```restart;
with(LinearAlgebra):
A := Matrix([[4, 2, -1], [2, 0, 2], [-1, 2, 0]], shape=symmetric, datatype=float):
x0:= <-1., 0., 1.>:
tol:=1e-4:
x1:=Normalize(A.x0, 2):
for i to 100 while Norm(x1-x0)>tol do
x0,x1 := x1, Normalize(A.x1, 2);
od:
x1^*.A.x1, x1, iter=i;
Eigenvalues(A); #check
```

## defective matrix...

Your matrix is defective (after plugging the conditions); for the 0 eigenvalue there are only 3 true eigenvectors.
You may want the Jordan form instead.

EIGENVECTOR-J.mw

## series...

A series is not designed to be used like this. It has a special structure, and you are messing with its operands.
Note that e.g.   2 + series(exp(x),x,3)  cannot be simplified by simplify

## Many solutions...

There are (too) many solutions. E.g.

```F:=(n,k) -> piecewise(k=2, n, k=n+1, n!, arbitrary);
```

F(F(n,n+1),2) assuming n<>1;
n!

It is also possible to write

```F:=proc(n,k)
if k=n+1 then return n! fi;
if k=2 then return n fi;
'F'(n,k) # arbitrary
end;```

In this case, assuming is not necessary.

## LinearAlgebra, Modular, Generic...

When the order of the field is prime, simply use LinearAlgebra:-Modular.
For a general finite field, use F:=GF(...) to define the field and then use the generic functions in the package LinearAlgebra:-Generic, e.g.  MatrixInverse[F](...).

## Just a bug...

The simplified version is:

```restart;
f:=k->Int(BesselJ(1, t)*BesselJ(0, k*t), t=0..infinity):
value(f(k));               # 1, wrong
value(f(k)) assuming k>1;  # 0, ok
value(f(k)) assuming k<1;  # 1, wrong (for k<=-1)
value(f(k)) assuming k<-1; # unevaluated,  why?
value([f(1/2),f(1),f(2)]) = evalf([f(1/2),f(1),f(2)]); # [1,1/2,0] , ok
```

You made me curious with this example.
(Actually it's obvious that you started from an approximation of ln(x)  [maybe using remez].)

So, you want the roots of f  in the interval (0,1]

It is possible to show that f has exactly 35 roots in 0..1 and they can be computed,
but unfortunately the Maple built-in commands are useless for this example.
Especially RootFinding:-NextZero  gave very poor results.
Disappointing!

```restart;
f := x -> 42312393/170170 + (87300630621*x)/5005 + (84260074354272*x^2)/1001 + (11572751542512000*x^3)/143 + (311217957498451500*x^4)/13 + (13736312974717541208*x^5)/5 + (702361109129611835904*x^6)/5 + (24407857262955082295808*x^7)/7 + (309128866743376578380625*x^8)/7 + (2034476230680567168673750*x^9)/7 + 971145727536841731347616*x^10 + (15981733578778631623709568*x^11)/11 + (3759286855176800298957060*x^12)/11 - (191464620989481690770115000*x^13)/143 - (1305974159036375354989560000*x^14)/1001 - (37785618862730180432031744*x^15)/91 - (8073200612643819138676179*x^16)/182 - (231388810612770205973655*x^17)/221 + (190600129650794094000*x^17 + 13377406242490734054600*x^16 + 195343515041861129011200*x^15 + 1040537189498550048000000*x^14 + 2518477612889131719000000*x^13 + 3093773724805440474252000*x^12 + 2048291569526360589849600*x^11 + 753720641133895782547200*x^10 + 155778902350755576750000*x^9 + 17974488732779489625000*x^8 + 1132669320760996761600*x^7 + 37459215415250044416*x^6 + 610748077422555072*x^5 + 4463801814780000*x^4 + 12619635840000*x^3 + 10816830720*x^2 + 1779084*x + 18)*ln(x):
plot(f, 1e-8.. 2*10^(-6))
```
```isolate(f(x), ln(x));
g:=unapply((lhs-rhs)(%), x);
```
```g1:=normal(diff(g(x),x)):
sturm(numer(g1),x,0,1),  sturm(denom(g1),x,0,1);
```

34, 0

# So, g' has 34 real roots and all are in (0,1);

```Digits:=100;
L:= [10.^(-6), fsolve(numer(g1)),  1.];
signum~(g~(L));nops(%);
# [-1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1,  -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, #  , -1, 1]
#                               36
# So, g (and f) has exectly 35 roots in (0,1);
```
`for i to 35 do i=fsolve(g, L[i]..L[i+1]) od;`

1 = .191561287749093307937227033160539875133752143080898731095086117213834851920515705439e-5
2 = .386863114821857362406592315527508687579950813861777293116197038157920498833807772632e-4
3 = .224368589492033177901033815331319924131620771048548665891126694787181114809420099728e-3
4 = .757628876465743496019898819891689451187953460530807054577709360768356026307038622925e-3
5 = .191586043560551126248587315094560734538794188822600052918532582221573212092426725234e-2
6 = .404341050752997444506457044613767563622770741756842016619777344000424267695214131236e-2
7 = .753803071564053290619849328637704614199820872264894616553942256990815977878634237764e-2
8 = .128347814103994790541831995624033323707358509644963097185620852545664106970618325968e-1
9 = .203881744484780444436159799720467561935237060222424328004944514734354971452828369636e-1
10 = .306531122499106815328715860139756718990911721033357327791034227500491192932066198881e-1
11 = .440652608872939466182326564099001428304047894621422939893923731352872694517742934307e-1
12 = .610215129442692726469217238624703315785764842073204122632174231202851842298289275127e-1
13 = .818611973531043298152537221149426167590700088952157881607430662062960493383888023674e-1
14 = 106848673577559840945862773565531573854056603799252847848710761297768054444573771543
15 = 136157907371746382188673002171784836243338456425810338155943739653948751569648355407
16 = 169859566116596553194242901269900419712969903139816354118873867296351077315658311599
17 = 207911095254057736286977980962266715752752350283085803618405408714970551082808317735
18 = 250150145986443098687961663490645213960126759514191210023876415880465673699179966025
19 = 296291621034620997624748934217305561232655042486558270046438747007382649214020765990
20 = 345928493082594386040296659177343956611809301967460342049504180924850028932218436257
21 = 398536433083337831429021148180533317041041457085873035259687560041254081758343858771
22 = 453482166534225841822476226178659839805335530836077276359472418786766336138346293288
23 = 510035358874306638935680974371652097595533372440279658507373385287879407012832470726
24 = 567383719969290157155679436652208052931514611589201379988400501404970528045712346664
25 = 624650915713400311490990374296520507166356685922856856979075726522784716531016773873
26 = 680916785283116652996592071240106370085097908173787711070599706365125792126746060999
27 = 735239288335241089539005060331848708632199246169167693175129656187672273735745985465
28 = 786677549795751681675503376421685045544117797715570060193731242723710977234911696913
29 = 834315332668459220663475591628563704363467776558796639367542621477772723238905506446
30 = 877284252847918207927215084652553361226192179200405713639782815744067137857243235030
31 = 914786055398803215378013441256970465033525746848992402044164144886322908035392340196
32 = 946113301720552053182303564465904800429669666884353849407725334749814613284598367425
33 = 970667885708695016778657128606495810300922043980656327177238682949045151388078173429
34 = 987977026592788760481681524535038678110293678929974123268967108570856367223442229329
35 = 997708817965433194699497427026838774265049425782232856894868033444326438291042212089

```r[0]:=1e-6;
for i do
r[i]:=RootFinding:-NextZero(g,r[i-1]+1e-10, maxdistance=1);
until r[i]=FAIL;      # disappointing!
```

# aborted !

Edit.  However, using
RootFinding:-NextZero(g,r[i-1]+1e-10, maxdistance=1, guardDigits=30);
NextZero works!

## Note that the parametrization is 4*Pi - ...

Note that the parametrization is 4*Pi - periodic in u, so, no need to plot over 0..6*Pi.

```restart;
f:=[arctan(-tan(u/4)) + ((cosh(v) - sinh(v))*sin(u))/2, -v/2 + ((cosh(v) - sinh(v))*cos(u))/2, 2*(cosh(v/2) - sinh(v/2))*sin(u/2)]:
e:=1.e-3:
P1:=plot3d(f, u = 0 .. 2*Pi-e,    v = -1 .. 1, grid = [80, 15], orientation = [50, 79], shading = XYZ, style = patch):
P2:=plot3d(f, u = 2*Pi+e .. 4*Pi, v = -1 .. 1, grid = [80, 15], orientation = [50, 79], shading = XYZ, style = patch):
plots:-display(P1,P2);
```