Jets.s 206 KB
Newer Older
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1
# J e t s 
Hynek Baran's avatar
Hynek Baran committed
2
#
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
3
# Differential calculus on jet spaces and diffieties
Hynek Baran's avatar
Hynek Baran committed
4
5
#

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
6
7
8
#
# See http://jets.math.slu.cz/ for further details
#
Hynek Baran's avatar
Hynek Baran committed
9

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
10
11
12
13
14
15
#
# A u t h o r s
#
# (c) 1999-2004 by Michal Marvan
# (c) 2005-2010 by Hynek Baran and Michal Marvan
#
Hynek Baran's avatar
Hynek Baran committed
16
17
18
19
20
21
22
23
#
#  L i c e n s e
#
# JETS is a free software distributed under the terms of the GNU General 
# Public License <http://www.gnu.org/copyleft/gpl.html> as published by 
# the Free Software Foundation.
#
# In particular, JETS comes with ABSOLUTELY NO WARRANTY.
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
24
25
26
27
28
#
# 
# A c k n o w l e d g m e n t s
#
# The first-named author was supported by GACR under project 201/07/P224
Hynek Baran's avatar
v 5.81    
Hynek Baran committed
29
# Developement of Jets was gratefully supported by project CZ.1.07/2.3.00/20.0002. 
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
30

Hynek Baran's avatar
Hynek Baran committed
31
32

#
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
33
34
# History of changes
#
Hynek Baran's avatar
v 5.71    
Hynek Baran committed
35
36
# v. 5.7  rel. 8 Dec 2010
# * Based on Jets code v. 5.6 as of 20 Jan 2010
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
37
#
Hynek Baran's avatar
v 5.71    
Hynek Baran committed
38
# v. 5.71 rel. 4 Jan 2011
Hynek Baran's avatar
v 5.72    
Hynek Baran committed
39
# * `divideout/unks/1`(): approved, better handling of polynomials
Hynek Baran's avatar
Hynek Baran committed
40
#
Hynek Baran's avatar
v 5.72    
Hynek Baran committed
41
42
43
44
# v 5.72
# * `index/traceless` matrix index function added
# * store(): formatting of output of assignments changed (each on separate line now)
# * new `jet/gen/all` and `count/gen/all` functions 
Hynek Baran's avatar
v 5.73a    
Hynek Baran committed
45
46
47
#
# v 5.73 alpha
# * first attempt to fix a critical bug in `cc/new` (some cc's was omitted by run)
Hynek Baran's avatar
Hynek Baran committed
48
49
50
51
52
53
54
55
56
57
#
# v 5.73 beta
# * `divideout/unks`, `type/divideout`: fixed (it was ok in very old jets but broken later)
#   and little improved ((exp(F)+1) is divided out now)
# * Varordering: new ordering revdeg implemented;
#   [degree, revdeg] is equivalent to the graded reverse lexicographic order (tdeg or grevlex)
# * `resolve/normalizer` added
# * Reporter implemented (quite primitive, resolve only)  
#
# v 5.74 alpha
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
58
59
# * resolve completly rewritten
# * started use of parallel MaP version of map (in resolve and derive) but no still not supported by Maple 16
Hynek Baran's avatar
Hynek Baran committed
60
61
62
# * derive redesigned
# * new functions Varl, VarL, LVar
# * minor evalTD bug fixed (wrong result on matrices) 
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
63
64
65
66
67
68
69
70
71
72
73
74
75
# * nonzero() bug (missing factor@) fixed (probably big issiue)
# * minor bug in Varordering fixed (if none fiber vars)
# * `resolve/nonresrat` limit introduced (price overlap of resolvable epr. over nonresolvable enforces FAIL of resolve)
# * minor improvements in reporting of assembled cc's
# * TestUnkOrd() introduced
# * sizefactor
# * `resolve/lin/price` back in the simplest form
# * `divideout/normalizer` added
# * Maple 14 compatibilty
#
#
# v 5.75 
# * old resolve is default
Hynek Baran's avatar
v 5.76    
Hynek Baran committed
76
77
#
# v. 5.76
Hynek Baran's avatar
v 5.77    
Hynek Baran committed
78
79
80
81
# * apply fixed (last ERROR replaced with unevaluated apply)
#
# v 5.77
# * `derive/1` : treating of polynoms and sums disabled (a bug found)
Hynek Baran's avatar
v 5.78    
Hynek Baran committed
82
83
84
85
86
87
88
#
# v 5.78
# package JetMachine introduced
# JetMachine:-Consequences introduced
# * `jet/gen/all` fixed
# * `convert/crack` introduced
# * reenabled Vars to accept sequential arguments
Hynek Baran's avatar
Hynek Baran committed
89
90
91
92
93
94
95
96
97
# * uncomment the line bellow to enable the new resolve implementation
# jets_new_resolve_enable := true;
#
#
# v 5.79
# * fixed missing local declaration in size() 
# * `resolve/nonresrat` reimplemented, resolve() FAILs when size condition is triggered
# * minor reporting improvements in `run/l`()
# * JetsProfiler introduced
Hynek Baran's avatar
v 5.81    
Hynek Baran committed
98

Hynek Baran's avatar
Hynek Baran committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
# * BasisExtractor function to symmetries/laws basis extraction introduced
# * `divideout/unks` improved (treating exponents)
# * `resolve/lin` unifyed
# * reporting in `resolve/lin` (lin resolvable/nonres. / nonlin)
# * CC:-markFF assertion commented out as a first aid, this issuie is to be resolved carefully
# * `TestUnkOrd/vars` weakened to report strong unknown ordering suspicients only
# * `run/l/extraders` introduced allowing to derive also subsets of cc in run (but none implementations besides empty)
# * fixed `AmICons/dependence` bug
# * fixed revdeg
# * putsize back
#
#
# v 5.80 beta
# * improved derive (noderives changed)
#
# 
# v 5.81
# * version mismatch cleanup
Hynek Baran's avatar
v 5.81    
Hynek Baran committed
117
118
119
120
121
122
# * apply -> vfapply M.M. 20 Feb. 2014
# * run() HB bugfix #if Bytes() > Blimit then reduce() fi
# * `divideout/unks`() bugfix: nonnegative numeric exponents only are divided out
# * jetorder() uses max()
# * `size/*/<`() at lest one result is always returned (low ressize and putsize settings do not cause error)
# * derive reporting slightly changed
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
123
124
125
#
#
# v 5.82
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
126
# * coordinates(), newfibre() and `jet/aliases`() has new option separator:="_" ()
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
127
# * using coordinates(..., separator="__"), jets are prettyprinted as subscripts (prior Maple 17)
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
128
# * `jet/aliases/mainseparator` is a global variant of the above option
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
129
# * TestUnkOrd() functionality limited to simple dependecy test (memory leak in second part of the test found)
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
130
# * BasisExtractor: unknown U collected in the result
131
132
133
134
#
#
# v 5.83
# * removed reduce() when store() is called (inside `store/pds`())
135
136
137
#
#
# v 5.84
Hynek Baran's avatar
Hynek Baran committed
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
# * fixed `resolve/lin` new implementation reportfail bug (jets_new_resolve_enable only)
# * linderive() introduced (but not used yet) 
#
#
# v 5.85
# * clear() rewritten:
#   + clear(pds) is also clearing cc's (unless `clear/pds/ccclearing`:='suppress' is assigned)
#   + clear(puts) and clear(assignments) introduced
#   + clear(..., output=eq) output format modifier
# * `puts/assignments`() and `put/name/tab` defined (as an generalisation of `unks/assignments`() resp. `unk/tab`)
# * store() is saving:
#   + all put symbols (form `put/name/tab`) instead of unknowns only (`unk/tab`)
#   + unknowns order (restored from `unk/tab`)
#   + varordering and Varordering
#   + terminal issue fixed  
#
#
# v 5.86
# * version mismatch (5.85) hopefully fixed 
# * `put/sizelimit` and `size/1/DD` introduced
# * run() refactored (cc's are derived)
# * `put/1` distinguishes redundand and contradictory attempt of putting already assigned expressions
Hynek Baran's avatar
v 5.87    
Hynek Baran committed
160
161
162
163
164
165
166
167
#
# v 5.87
# * not deriving cc's by default (`run/l/extraders` used in slightly different mode)
# * `run/1` reporting cleanup
# * Report Reportf macros rewritten, rt, rb no more needed
# * `size/1/DD` abandoned for a hidden bug
# * `put/limit/length`, `put/limit/size` instead `put/sizelimit` 
# * JetMachine[Consequences]:-AmICons: Using `AmICons/ignore`() for nonzero(), Varordering() and unknowns()
Hynek Baran's avatar
v 5.89    
Hynek Baran committed
168
#
Hynek Baran's avatar
v 5.88    
Hynek Baran committed
169
170
171
172
173
# v 5.88
# * fixed `resolve/lin` new implementation reportfail bug (jets_new_resolve_enable only)
# * linderive() introduced (but not used yet) 
# * `unks/TD` uses forceError=true in `vars/1` calls
# * JetMachine[Consequences] testing script changed (currentdir)
Hynek Baran's avatar
v 5.89    
Hynek Baran committed
174
175
# * `clear/assignments` is clearing `put/name/tab` properly# current
# * JetMachine[Consequences] testing script changed (currentdir)
Hynek Baran's avatar
v 5.88    
Hynek Baran committed
176
# * `clear/assignments` is clearing `put/name/tab` properly
Hynek Baran's avatar
v 5.89    
Hynek Baran committed
177
178
179
#
# v 5.89
# * run: linear eqs are ALL passed to resolve
Hynek Baran's avatar
Hynek Baran committed
180
181
182
183
#
# v 5.90
# * jet: option remember bug (multiple eqs problem as reported by JSK) fixed by introducing `eqn/table`
# * lengthselect introduced
Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
184
185
186
187
188
189
190
#
# v. 5.91
# * Vars has option cache`
# * resolve FAILs collecting introduced
# * NewIntSeq() introduced
# * Report macro changed to inlined function
# * New resolve implementation (jets_new_resolve_enable) removed from the Jets.s source file  to the Jets.newresolve.s file
Hynek Baran's avatar
Hynek Baran committed
191
#
Hynek Baran's avatar
v 5.93    
Hynek Baran committed
192
# v 5.93
Hynek Baran's avatar
Hynek Baran committed
193
# * testing new nonlinear resolve implementation in file (cat(jets_new_resolve_enable, "/Jets.newresolve.s"))
194
# * LVar returns NULL on expressions without Vars
195
# * numerous bug fixes in the computation of coverings (MM 23.4.2019)
Hynek Baran's avatar
v 5.88    
Hynek Baran committed
196

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
197
198
199
200
201
202
203
###########################################################################################
###########################################################################################
# Initialization and configuration
###########################################################################################
###########################################################################################

interface(screenwidth=120):
Hynek Baran's avatar
v 5.93    
Hynek Baran committed
204
printf("Jets 5.93  as of May 2, 2019 (%s)\n",
Hynek Baran's avatar
Hynek Baran committed
205
      StringTools[FormatTime]("%Y-%m-%d %T", timestamp=FileTools[Status](__FILE__)[5]));
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
206
207
208
209
210
211
212
213
214
215
216

#
# Source code configuration, options and parameters
#

### optional features
# Some parts of Jets source code is disbbled by default,
# since it is not proprely finished or tested yet.
# Use assignments bellow (before reading Jets.s in) to enable such additional code.

#`Jets/opts`["Optionals"]["Multiord"] := true:
Hynek Baran's avatar
Hynek Baran committed
217
#jets_new_resolve_enable := ".";
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265

#
# Debugging
#

### reporting
#
# Reporting about the progression of calculation is available independently of debug mode bellow.
# Syntax: reporting(derive = 0, resolve = 0, run = 0, cc = 0, pd = 0); 

### debug mode switch
#
# Some debugging functionality (as fuction arguments type checks, assumptions or logging) 
# is disabled by default but may be turned on by the macro bellow. 
# A lot of computing time is saved when debug mode is disabled, 
# but the AssertLevel() function and some more debugging functionality is not available.
# Turning debugging mode on may dramatically harm the time and memory usage.
#
# In the case of any problem (as an error raised or wrong result returned),
# please turn the debugging mode ON by removing the very first character # in the line bellow:
#
#$define jets_mode_debugging # <<<- TURN ON DEBUG MODE HERE by UNcommenting this line #  <- DEBUG ON/OFF


### debugging code
$ifdef jets_mode_debugging
# the case when debugging is ON:
print(`debugging mode enabled`);
print(kernelopts(version));
kernelopts(opaquemodules=false): 
$define ATC(a,t) a::t # Arguments Type Control enabled
$define DOE(a) a      # Debug Only Expressions enabled
$else
# the case when debugging is OFF:
$define ATC(a,t) a # Arguments Type Control disabled
$define DOE(a)     # Debug Only Expressions ignored
$endif

#
# Initialization tests
#

### single load test
# make sure that this source file is not already loaded:
if assigned(cat(`jets/read/flag/`,__FILE__)) then error "source file %1 already loaded.", __FILE__; fi;
if assigned(cat(`jets/read/flag/`,__FILE__)) then `quit`(255); fi; # force immediate exit from read() statement
assign(cat(`jets/read/flag/`,__FILE__), true):

Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
266
### Reporting macros
267
268
269
270
271

ProcNameN := proc(n) option inline; debugopts('callstack')[1+1+3*n] end: # The list will contain 3n+4 elements. The first element is the name DEBUGSTACK. The remaining elements are in groups of 3.
ProcBaseNameN  := proc(n) option inline; StringTools:-Split(convert((ProcNameN(n)), string), "/")[1] end:
#ProcBaseName   := proc() option inline; remove(`=`, map(ProcBaseNameN, [$ 1..((nops(debugopts('callstack'))-1)/3)]), "unknown")[1] end: ### VERY slow
ProcBaseName   := proc() option inline; ProcBaseNameN(0) end:
Hynek Baran's avatar
v 5.87    
Hynek Baran committed
272
273
ProcBaseSymbol := proc() option inline; convert(ProcBaseName(), symbol) end:
ProcBaseSYMBOL := proc() option inline; convert(StringTools:-UpperCase(ProcBaseName()), symbol) end:
Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
274
275
276
277

Report := proc(l, m)
  option inline;   
  `if` (`report/tab`[ProcBaseSymbol()] > l, # l<=0 or ... for some magic reason wont work
Hynek Baran's avatar
Hynek Baran committed
278
        report(cat(ProcBaseSYMBOL(), '`:`'), [cat(ProcNameN(0),'`[`',l,'`]`',":"), `if`(type(m,list),op(m),m)]),
Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
279
280
281
282
283
284
        NULL);
end:

Reportf := proc(l, m)
  option inline;  
  `if`(`report/tab`[ProcBaseSymbol()] > l,
Hynek Baran's avatar
Hynek Baran committed
285
       report(cat(ProcBaseSYMBOL(), '`:`'), sprintf(cat("%a[%a]: ", op(1,m)), ProcNameN(0), l, op(2..-1, m))),
Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
286
287
288
289
       NULL);
end:


Hynek Baran's avatar
Hynek Baran committed
290
291
292

##############################
#MaP := evalf(Threads[Map]);
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
293
MaP := evalf(map):
Hynek Baran's avatar
Hynek Baran committed
294

Hynek Baran's avatar
v 5.75    
Hynek Baran committed
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
#
# old Maple version backward compatibility
#
# the routines defined in newer versions only are called as compat[routine]
# in the code bellow are these routines 
# - assigned to this auxiliary table (if available in the current version)
# - implemented by Hynek for backward compatibility 

### Maple 14 is probably the minimum
if not(assigned(DifferentialAlgebra)) then # DifferentialAlgebra available since Maple 14
  lprint(kernelopts(version));
  error "Maple older than 14 detected, Jets cannot be used.";
fi;
if not(assigned(DifferentialAlgebra)) then
  `quit`(255);
fi;
### Maple 15 features:
if not(assigned(numelems)) then # numelems available since Maple 15
  printf("Warning: Maple of version older than 15 detected. Please upgrade to Maple 15.\n%s\n", kernelopts(version));
  compat[Record[packed]]:= Record: # Record[packed] not available, use Record instead
else
  # we have Maple 15 (or newer)
  compat[Record[packed]] := Record[packed]:
fi:
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
319

Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
#
# Aux macros
#

# create generator of integer sequence
NewIntSeq := proc(initvalue::integer:=0) 
  local i := initvalue;
  proc({`set`::{integer,identical(NULL)}:=NULL, `get`::truefalse:=false})
    if type(`set`, integer) then i := `set`;
    elif `get` then return i;
    else i := i+1
    fi;
  end
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
335
336
337
338
339
340
341
###########################################################################################
###########################################################################################
# Jets - Main procedures
###########################################################################################
###########################################################################################

# This is a set of main routines based on the original Jets.s by MM
Hynek Baran's avatar
Hynek Baran committed
342
343
344
345
346
347
348
349
350
351
352
353
354

#
#  D e c l a r a t i o n s
#

# A procedure whose purpose is to change settings of global variables 
# is called a declaration.
# If named in plural, it resets the previous declaration.
# If named in singular, it has a cummulative effect.

#
#  A u x i l i a r y   p r o c e d u r e s
#
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
355
356
357

Existing := proc(c,a) assigned(c||a) end:
Call := proc(c,a) c||a end: # || replaces . since 6.0
Hynek Baran's avatar
Hynek Baran committed
358
    
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
359
360
# sorting by size

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
361
size := proc() # HB generalized size of any number of arguments
Hynek Baran's avatar
Hynek Baran committed
362
  local a;
Hynek Baran's avatar
Hynek Baran committed
363
  #try 
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
364
    add(`size/1`(a), a in [args]);
Hynek Baran's avatar
Hynek Baran committed
365
366
367
368
369
  #catch:
  #  printf("\nsize: failed and returned infinity: %q\n", 
  #         StringTools[FormatMessage](lastexception[2..-1]));
  #  infinity;
  #end:
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
370
371
end:

Hynek Baran's avatar
v 5.75    
Hynek Baran committed
372
sizefactor := 1000000000:
Hynek Baran's avatar
Hynek Baran committed
373

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
374
`size/1` := proc(a) # MM size evaluator
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
375
  global sizefactor;
Hynek Baran's avatar
Hynek Baran committed
376
  local b;
Hynek Baran's avatar
v 5.87    
Hynek Baran committed
377
378
  #b := eval(eval(a,  TD=`size/1/DD`), pd=`size/1/DD`);
  evalf(nops(Vars(a)) + length(a)/sizefactor);
Hynek Baran's avatar
Hynek Baran committed
379
380
end:

Hynek Baran's avatar
v 5.87    
Hynek Baran committed
381
382
383
384
385
386
387
388
## there was a hidden bug, consider
## 1/(pd(F,u)-pd(F,v)) --> 1/(DD(F,d)-DD(F,d)) = 1/0
##
#`size/1/DD` := proc(f, p)
#  # do not count too much of TD(f,p), pd(f,p) derivating monomial p
#	local d; # dummy 
#	'DD'(f, cat(d $ degree(p))) # slight penalty for higher derivatives
#end:
Hynek Baran's avatar
Hynek Baran committed
389

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
390
391
392
393
394
395
396
#sizesort := proc(ql,pr)  # ql = list of expressions, pr = sizeing proc
#  local qls;
#  qls := map(proc(q,pr) [q,pr(q)] end, ql, pr);
#  qls := sort(qls, proc(a,b) evalb(op(2,a) < op(2,b)) end);
#  map (proc(a) op(1,a) end, qls)
#end:
sizesort := proc(L::list, sizef)
Hynek Baran's avatar
v 5.87    
Hynek Baran committed
397
398
  local l;
  return map(attributes,sort([seq](setattribute(SFloat(sizef(l),0),l), l=L),`<`));
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
399
400
end proc:
# see http://www.mapleprimes.com/blog/joe-riel/sorting-with-attributes
Hynek Baran's avatar
Hynek Baran committed
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420

# reducing products

reduceprod :=  proc(a)   
  if type (a,`^`) then
    if type (op(2,a), posint) then procname(op(1,a))
    elif type (op(2,a), negint) then 1
    else a
    fi
  elif type (a,`*`) then convert(map(procname,[op(a)]),`*`)
  else a
  fi
end:

#
#  V a r i a b l e s 
#

# Declaration of coordinates (only plural form):

Hynek Baran's avatar
v 5.82    
Hynek Baran committed
421
coordinates := proc (blist,flist, MaxAliasDegree::integer:=0, {separator::string:=`jet/aliases/mainseparator`})
Hynek Baran's avatar
Hynek Baran committed
422
423
  global `b/var/list`,`b/var/s`,`b/dim`,`b/<</list`,
    `f/var/list`,`f/var/s`,`f/dim`,`f/<</list`;
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
424
  local res;
Hynek Baran's avatar
Hynek Baran committed
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
  if nargs = 0 then
    ERROR(`arguments should be:\n
[list of base variables], [list of fibre variables], optional maxorder`) 
  fi;
  if not type([op(blist)], list(name)) then
    ERROR(`base coordinates must be unassigned names`)
  fi;
  if not type([op(flist)], list(name)) then
    ERROR(`fibre coordinates must be unassigned names`)
  fi;
  if nops([op(blist),op(flist)]) <> nops({op(blist),op(flist)}) then
    ERROR(`coordinates must be mutually different`)
  fi;
  `b/var/list` := blist;
  `b/var/s` := {op(blist)};
  `b/dim` := nops(`b/var/s`);
  `b/<</list` := `b/var/list`:
  `f/var/list` := flist;
  `f/var/s` := {op(flist)};
  `f/dim` := nops(`f/var/s`);
  `f/<</list` := `f/var/list`;
  noderives():
Hynek Baran's avatar
Hynek Baran committed
447
  doderives();
Hynek Baran's avatar
Hynek Baran committed
448
449
  refresh ();
# Result:
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
450
451
452
453
454
455
  if MaxAliasDegree > 0 then 
    res := `jet/aliases`(`f/var/list`, MaxAliasDegree, ':-separator'=separator);
    if separator <> "_" then # for backward compatibility, allways alias u_x
       `jet/aliases`(`f/var/list`, MaxAliasDegree, ':-separator'="_") 
    fi; 
    res;
Hynek Baran's avatar
Hynek Baran committed
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
  else `b/var/list`, `f/var/list`
  fi
end:

`b/var/list` := []: `b/var/s` := {}: `b/dim` := 0:
`f/var/list` := []: `f/var/s` := {}: `f/dim` := 0:

# Type checking procedures:

`type/b/var` := proc (x) member (x, `b/var/s`) end:
`type/f/var` := proc (f) member (f, `f/var/s`) end:

# Jet variables are unassigned calls to 'jet'.

`type/j/var` := specfunc(anything,jet):

# Variables are of type 'var'.

`type/var` := proc(x)
  if type (x,'name') then type (x,{`b/var`, `f/var`})
  else type (x,specfunc(anything,jet))
  fi
end:

#
#  C o u n t s 
#

# A count is a product of nonnegative integer powers, but 1 is not a count.

`type/count` := proc (x) 
  local x1;  
  if type (x,`^`) then type (op(2,x), posint) 
  elif type (x,`*`) then
    for x1 in x do if not type (x1,'count') then RETURN(false) fi od;
    true
  elif x = 1 then false
  else true
  fi
end:

# A count is proper when it is a product or a power

`type/count/proper` := {`*`,`^`}:

# A count is of first order when it is neither a product nor a power

`type/count/1order` := proc(x) not type(x,`count/proper`) end:

# A b/count is a product of nonnegative integer powers of base variables,
# but 1 is not a b/count.

`type/b/count` := proc (x) 
  local x1;  
  if type (x,`^`) then type (op(1,x),`b/var`) and type (op(2,x), posint) 
  elif type (x,`*`) then
    for x1 in x do if not type (x1,`b/count`) then RETURN(false) fi od;
    true
  else type (x,`b/var`)
  fi
end:

# Dealing with counts: 
# `count/f`(x) := the first ``factor'' of the count x;
# `count/r`(x) := x/`count/f`(x);
# `count/length`(x) := the sum of the powers in x;
# `count/sgn`(x) := the parity of `count/length`x;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
523
# `transform/count`(x) :=  the product of non-negative powers in x # = numerator 
Hynek Baran's avatar
Hynek Baran committed
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
# (transforms x into 1 or count).

`count/f` := proc (x)
  if type (x,`^`) then op(1,x)
  elif type (x,`*`) then `count/f`(op(1,x))
  else x
  fi
end:

`count/r` := proc (x)
  if type (x,`^`) then x/op(1,x)
  elif type (x,`*`) then x/`count/f`(op(1,x))
  else NULL
  fi
end:

`count/length` := proc (x)
  if x = 1 then 0 
  elif type (x,`^`) then op(2,x)
  elif type (x,`*`) then convert (map (`count/length`, [op(x)]), `+`)
  else 1
  fi
end:

`count/sgn` := proc (x) 1 - 2*(`count/length`(x) mod 2) end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
550
551
552
553
554
555
556
557
558
#`transform/count` := proc(x)
#  local aux;
#  if type (x,`^`) then if op(2,x) < 0 then 1 else x fi
#  elif type (x,`*`) then aux := map(`transform/count`, x);
#  else x
#  fi
#end:

`transform/count` := op(numer):
Hynek Baran's avatar
Hynek Baran committed
559

Hynek Baran's avatar
v 5.73a    
Hynek Baran committed
560
`count/gen/all` := proc(b::list, m::{posint,0}) # HB
Hynek Baran's avatar
v 5.72    
Hynek Baran committed
561
562
563
564
565
566
567
568
569
570
571
572
573
  description "Generates all counts in variables b up to m-th order";
  `count/gen/all/1`(b,m)[2..-1]; # 1 is not a count
end:

`count/gen/all/1` := proc(b::list, m::{posint,0}) # HB
  if nops(b) = 0 or m<=0 then 
    1
  else
    seq(seq(i*j, j in [procname(b[2..-1], m-`count/length`(i))]), i in seq(b[1]^i, i=0..m));
  fi
end:


Hynek Baran's avatar
Hynek Baran committed
574
575
576
577
578
579
580
581
582
# Jet order of a variable q

varorder := proc(q)
  if type(q,'name') then 0
  else `count/length`(op(2,q))
  fi
end:

jetorder := proc(f)
Hynek Baran's avatar
v 5.81    
Hynek Baran committed
583
  max(jetorders(f)) # HB
Hynek Baran's avatar
Hynek Baran committed
584
585
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
586
jetorders := proc(f) # HB
Hynek Baran's avatar
v 5.81    
Hynek Baran committed
587
  map(varorder,`vars/1`(f))
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
588
589
end:

Hynek Baran's avatar
v 5.82    
Hynek Baran committed
590
591
592
593
594
595
596
597
598
599
600
601
602
# If you assign
# `jet/aliases/mainseparator` := "__" ;
# BEFORE reading-in jets.s
# jets will be pretyprinted using subscripts, i. e. aliased using "__" instead of "_"
# jet(u,x) = u__x
# and old fashion jet aliases will be still available for backward compatibility
# jet(u,x)=u__x=u_x
if not(assigned(`jet/aliases/mainseparator`)) then `jet/aliases/mainseparator` := "_" fi:

# Creating aliases for jet variables. Arguments are:
# f1,...,fn,r, where r = jet order and 
# f's are fibre variables (optional)

Hynek Baran's avatar
Hynek Baran committed
603
604
605
606
# Creating aliases for jet variables. Arguments are:
# f1,...,fn,r, where r = jet order and 
# f's are fibre variables (optional)

Hynek Baran's avatar
v 5.82    
Hynek Baran committed
607
`jet/aliases` := proc ({separator::string:="_"})
Hynek Baran's avatar
Hynek Baran committed
608
609
  global `b/var/list`;
  local aux,flist,r;
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
610
611
612
  r := _rest[-1];
  if _nrest = 1 then flist := `f/var/list` 
  else flist := _rest[1..-2] 
Hynek Baran's avatar
Hynek Baran committed
613
614
615
616
  fi;
  if not type (`b/var/list`,'list'(symbol)) then
    ERROR (`base variables must be symbols`) 
  fi;
Hynek Baran's avatar
v 5.72    
Hynek Baran committed
617
618
619
  #aux := expand((1 + convert(`b/var/s`,`+`))^r - 1);  
  #aux := map(proc(a) a/coeffs(a) end, [op(sort(aux, `b/var/list`))]);
  aux := sort([`count/gen/all`(`b/var/list`, r)]); # HB
Hynek Baran's avatar
Hynek Baran committed
620
  alias(op(map(proc(u,aux) if type(u,symbol) then
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
621
        op(map(proc(x,u) cat(u, separator, `bcount//str`(x, separator)) = 'jet'(u,x) end,
Hynek Baran's avatar
Hynek Baran committed
622
623
624
          aux, u))
      else
        op(map(proc(x,u)
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
625
            cat(op(0,u), separator, `bcount//str`(x, separator))[op(u)] = 'jet'(u,x)
Hynek Baran's avatar
Hynek Baran committed
626
627
628
629
630
          end, aux, u))
      fi
    end, flist, aux)));
end:

Hynek Baran's avatar
v 5.82    
Hynek Baran committed
631
`bcount//str` := proc (x, separator)
Hynek Baran's avatar
Hynek Baran committed
632
633
634
  if type (x,'name') then x
  elif type (x,`^`) then 
    if op(2,x) > 3 then cat(op(2,x),op(1,x)) else op(1,x) $ op(2,x) fi
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
635
  elif type (x,`*`) then `bcount//str/_`(op(map(`bcount//str`, [op(x)], separator)), ':-separator'=separator)
Hynek Baran's avatar
Hynek Baran committed
636
637
638
639
  else ERROR (`not a count`, x)
  fi 
end:

Hynek Baran's avatar
v 5.82    
Hynek Baran committed
640
`bcount//str/_` := proc ({separator::string:="_"})
Hynek Baran's avatar
Hynek Baran committed
641
  local i,a,b,ans;
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
642
643
644
645
  ans := _rest[1];
  for i from 2 to _nrest do
    a := _rest[i-1]; b := _rest[i];
    if length(a) + length(b) > 2 then ans := ans,separator,b else ans := ans,b fi
Hynek Baran's avatar
Hynek Baran committed
646
647
648
649
650
651
652
  od 
end:

`count//seq` := proc (x)
  if type (x,`^`) then `count//seq`(op(1,x)) $ op(2,x)
  elif type (x,`*`) then op(map(`count//seq`,[op(x)]))
#  elif type (x,symbol) then x
Hynek Baran's avatar
v 5.82    
Hynek Baran committed
653
#  else cat(op(1,x),'_',`bcount//str`(op(2,x), separator)) 
Hynek Baran's avatar
Hynek Baran committed
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
  else x
  fi
end:

#
#  N a m e s
#

`name/tab` := table([]):

# Unassigned names may acquire a meaning through a declaration.
# Such names are said to be registered.
# For example, variables are registered names since they are declared
# through the coordinates command.

# Meanings other than 'variable' are stored in `name/tab`.
# registered() prints meanings of all names;
# registered(a) prints all names of meaning a;

# WARNING: After  a s s u m e,  names lose their registration and must be
# registered again.

registered := proc() 
  local a;
  if nargs = 0 then
    op(select(proc(e) op(1,e) = eval(op(1,e)) end, op(2,op(`name/tab`))))
  else a := args[1];
    op(map(proc(e,a) op(1,e) end, select(
      proc(e,a)
        op(1,e) = eval(op(1,e)) and type(op(1,e),a) 
      end, op(2,op(`name/tab`)), a), a))
  fi
end:

# To clear a declaration from all names:

Hynek Baran's avatar
Hynek Baran committed
690
691
clear := proc(opts::seq(name), {output::identical(expr, eq):=expr})
  local a,aux,ans,nex;
Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
692
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
  aux := {opts};
  #lprint(aux);
  #if aux minus {} <> {args} then
  #  ERROR(`not a name`, op({args} minus aux)) 
  #fi;
  
  if has(aux, puts) then aux := aux union {cc, pds, assignments}
  elif has(aux, pds) and `clear/pds/ccclearing`<>'suppress' then aux := aux union {cc} fi;
  
  aux, nex := selectremove(proc(a) Existing(`clear/`,a) end, aux);
  
  #lprint(aux,nex);
  #if aux <> {args} then
  #  ERROR(`invalid argument`, op({args} minus aux)) 
  #fi;
Hynek Baran's avatar
Hynek Baran committed
708
  ans := NULL;
Hynek Baran's avatar
Hynek Baran committed
709
710
711
  aux := sort(convert(aux,list)); # assignments is the first in the alphabet :)
  for a in aux do ans := ans, Call(`clear/`,a)(_options) od;
  return [ans];
Hynek Baran's avatar
Hynek Baran committed
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
end:

#
#   P a r a m e t e r s
#

parameter := proc()
  global `name/tab`;
  local a,aux;
  aux := select (type,{args},name);
  if aux <> {args} then ERROR(`not a name`, op({args} minus aux)) fi;
  aux := select(
    proc(a)
      type (a,'dep')
      or type (a,{`b/var`,`f/var`,'constant'})
      or (type (`name/tab`[a],'set')
        and `name/tab`[a] minus {'parameter','nonlocal'} <> {})
    end, {args});
  if aux <> {} then ERROR (`name(s) already used`, op(aux)) fi;
  for a in [args] do 
    if type (`name/tab`[a],'set') then
      `name/tab`[a] := `name/tab`[a] union {'parameter'}
    else `name/tab`[a] := {'parameter'}
    fi
  od;
  registered('parameter')
end:

`clear/parameter` := proc()
  local a,aux;
  global `name/tab`;
Hynek Baran's avatar
Hynek Baran committed
743
  print(`Clearing parameter's.`);
Hynek Baran's avatar
Hynek Baran committed
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
  aux := {registered('parameter')};
  for a in aux do
    `name/tab`[a] := `name/tab`[a] minus {'parameter'} 
  od;
  op(aux)
end:

parameters := proc()
  `clear/parameter`();
  parameter(args)
end:

# Parameters are of type 'parameter':

`type/parameter` := proc (x)
  if type (`name/tab`[x],'set') then
    member ('parameter',`name/tab`[x])
  else false
  fi
end: 

# Type 'ar' is either 'parameter' or 'var':

`type/ar` := {'var','parameter'}:

`type/ar/count` := proc (x) 
  local x1;  
  if x = 1 then false 
  elif type (x,`^`) then type (op(1,x),'ar') and type (op(2,x), posint) 
  elif type (x,`*`) then
    for x1 in x do if not type (x1,`ar/count`) then RETURN(false) fi od;
    true
  else type (x,'ar')
  fi
end:

#
#   D e p e n d e n c e
#

# Names may depend on variables:

dependence := proc ()
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
787
  local f,fs,es,e, opts; 
Hynek Baran's avatar
Hynek Baran committed
788
  global `dep/tab`;
Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
789
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
790
791
  fs := select(type, {args}, function);
  es := select(type, {args}, `=`);
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
792
793
  opts := select(type, {args}, symbol); # HB
  if fs union es union opts <> {args} then # HB
Hynek Baran's avatar
Hynek Baran committed
794
795
796
797
798
799
800
801
802
803
    ERROR (`wrong arguments`, op({args} minus fs minus es)) 
  fi; 
  es := es union map(proc(f) op(0,f) = {op(f)} end, fs);
  fs := select(
    proc(e)
      type (op(1,e),{`b/var`,`f/var`,'constant','parameter'})
    end, es);
  if fs <> {} then
    ERROR (`name(s) already used`, op(fs))
  fi; 
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
804

Hynek Baran's avatar
Hynek Baran committed
805
806
807
808
809
810
811
812
  for e in es do 
    if type (op(2,e), set('ar')) then `dep/tab`[op(1,e)] := op(2,e)
    else ERROR (`forbidden type of dependence`)
    fi
  od;
  refresh(); 
  es := {op(op(2,op(`dep/tab`)))};
  es := select(proc(e) evalb(op(1,e) = eval(op(1,e))) end, es);
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
813
814
815
816
817
818
819
  # HB:
  if `subs` in opts then
    seq((op(1,e)=op(1,e)(op(op(2,e)))), e = es) 
  else  
  # :HB
    seq(op(1,e) = op(2,e), e = es)
  fi   
Hynek Baran's avatar
Hynek Baran committed
820
821
822
823
824
825
826
827
828
end:

`dep/tab` := table ():

`type/dep` := proc(a) type(`dep/tab`[a],'set') end:

`clear/dependence` := proc()
  local e,el,ans;
  global `dep/tab`;
Hynek Baran's avatar
Hynek Baran committed
829
  print(`Clearing dependence's.`);
Hynek Baran's avatar
Hynek Baran committed
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
  el := op(2,op(`dep/tab`));
  ans := {seq((op(1,e))(op(op(2,e))), e = el)};
  `dep/tab` := table([]);
  ans 
end:

dependences := proc()
  `clear/dependence`();
  dependence(args)
end:

#
#   T o t a l   d e r i v a t i v e
#

TD := proc (f, x::`b/count`)
  global `EVALTD/s`,`EVALTD/b`,`EVALTD/n`;
  if nargs = 1 then RETURN (`TD~`(procname,f)) fi;
  if nargs > 2 then RETURN (TD(TD(f, args[2..nargs-1]), args[nargs])) fi;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
849
  if type (f,{'numeric', 'complex'}) then 0 # HB 5. 5. 2008
Hynek Baran's avatar
Hynek Baran committed
850
851
852
853
854
855
  elif type (f,'name') then
    if member (f,[constants]) then 0
    elif type (f,`b/var`) then if f = x then 1 else 0 fi
    elif type (f,`f/var`) then jet(f,x)
    elif type (f,'parameter') then 0
    elif type (f,'dep') then
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
856
      if `vars/1`(f) minus `b/var/s` = {} then pd(f,x)
Hynek Baran's avatar
Hynek Baran committed
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
      elif not `EVALTD/b` and not member(f,`EVALTD/s`)
        and `count/length`(x) > `EVALTD/n` then 'TD'(f,x) 
      elif type (x,'name') then 
        if type (f,'dep') then `TD/dep/s`(f,x,`dep/tab`[f]) fi
      else TD (TD(f,`count/f`(x)), `count/r`(x)) 
      fi
    elif Existing(`TD/indexed/`,op(0,f)) then
      Call(`TD/indexed/`,op(0,f))(op(f),x)
    else 'TD'(f,x)
    fi
  elif type (f,`+`) then map (procname, f, x)  
  elif type (f,'function') then 
    if type (f, `j/var`) then jet(op(1,f),x*op(2,f)) 
    elif type (f, specfunc(anything,pd)) then
      if type (op(1,f),'dep') then 
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
872
        if `vars/1`(op(1,f)) minus `b/var/s` = {} then pd(op(1,f), x*op(2,f))
Hynek Baran's avatar
Hynek Baran committed
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
        elif not `EVALTD/b` and not member (op(1,f),`EVALTD/s`)
          and `count/length`(x) > `EVALTD/n` then 'TD'(f,x)
        elif type (x,'name') then `TD/dep/s`(f,x,`dep/tab`[op(1,f)]) 
        else TD (TD(f,`count/f`(x)), `count/r`(x)) 
        fi
      else 'TD'(f,x)
      fi
    elif type (f, specfunc(anything,TD)) then 'TD'(op(1,f), x*op(2,f))
    elif type (x,'name') then
      if Existing(`der/`,op(0,f)) then
        Call(`der/`,op(0,f))(procname,op(f),x)
      else ERROR (`unexpected function`, f) 
      fi
    else TD (TD(f,`count/f`(x)), `count/r`(x))
    fi
  elif type (x,'name') then
    if type (f,`*`) then `der/*` (procname, f, x)
    elif type (f,`^`) then `der/^` (procname, op(1,f), op(2,f), x)
    else ERROR (`unexpected object`, f) 
    fi
  else TD (TD(f,`count/f`(x)), `count/r`(x))
  fi
end:

`TD~` := proc (pn,f)
  if type (pn,'indexed') then 
    if nops(pn) <> 1 then ERROR (`too many indices`, op(pn)) 
    else RETURN (TD(f,op(1,pn)))
    fi
  else eval(f)
  fi
end:

`TD/dep/s` := proc (f,x,s)
  convert(map(proc(p,f,x) pd(f,p)*TD(p,x) end, [op(s)], f,x), `+`)
end:

# Evaluating total derivatives

evalTD := proc (f) 
  local a,b,`EVALTD/b/old`,`EVALTD/s/old`;
  global `EVALTD/b`, `EVALTD/s`;
  `EVALTD/b/old` := `EVALTD/b`; `EVALTD/s/old` := `EVALTD/s`;
Hynek Baran's avatar
Hynek Baran committed
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
  if type(f,sequential) or type(f, matrix) then # HB
    map('procname', f, args[2..-1])
  else
    if nargs = 1 then `EVALTD/b` := true
		elif type (args[2], set(name)) then `EVALTD/b` := false;
		 `EVALTD/s` := args[2]
		else ERROR (`second argument must be a set of names`)
		fi; 
		b := eval(f); 
		a := {}; # a either different from a or evalTD(a) = a
		while hasfun(b,TD) and b <> a do
			a := b;
			b := traperror(eval(a))
		od;
		`EVALTD/b` := `EVALTD/b/old`; `EVALTD/s` := `EVALTD/s/old`;
		if b = lasterror then ERROR (lasterror) fi;
		b
  fi;
Hynek Baran's avatar
Hynek Baran committed
934
935
end: 

Hynek Baran's avatar
Hynek Baran committed
936
937


Hynek Baran's avatar
Hynek Baran committed
938
939
940
941
942
943
944
945
946
947
948
`EVALTD/s` := {}:
`EVALTD/b` := false:
`EVALTD/n` := 0:

#
#   P a r t i a l   d e r i v a t i v e s
#

pd := proc (f, p::`ar/count`)
  if nargs = 1 then RETURN (`pd~`(procname,f)) fi;
  if nargs > 2 then RETURN(pd(pd(f, args[2..nargs-1]), args[nargs])) fi;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
949
  if type (f,{'numeric', 'complex'}) then 0 # HB 5. 5. 2008
Hynek Baran's avatar
Hynek Baran committed
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
  elif type (f,'name') then
    if member (f,[constants]) then 0
    elif type (f,'ar') then if f = p then 1 else 0 fi
    elif Existing (`pd/indexed/`,op(0,f)) then
      Call(`pd/indexed/`,op(0,f))(op(f),p)
    else `pd//tab`(f,p)
    fi
  elif type (f,`+`) then map (procname, f, p)  
  elif type (f,'function') then 
    if type (f, `j/var`) then if f = p then 1 else 0 fi
    elif type (f, specfunc(anything,pd)) then `pd//tab`(op(1,f), op(2,f)*p) 
    elif type (p,'ar') then
      if type (f, specfunc(anything,TD)) then `pd/TD`(op(f), p)     
      elif Existing(`der/`,op(0,f)) then
        Call(`der/`,op(0,f)) (procname,op(f),p)
      else `der//unknown`(procname,op(0,f),op(f),p)  # M.M. 26.11.2004
      fi
    else pd (pd(f,`count/f`(p)), `count/r`(p))
    fi
  elif type (p,'ar') then
    if type (f,`*`) then `der/*` (procname, f, p)
    elif type (f,`^`) then `der/^` (procname, op(1,f), op(2,f), p)
    else ERROR (`unexpected object`, f) 
    fi
  else pd (pd(f,`count/f`(p)), `count/r`(p))
  fi
end:

`pd~` := proc (pn,f)
  if type (pn,'indexed') then 
    if nops(pn) <> 1 then ERROR (`too many indices`, op(pn)) 
    else RETURN (pd(f,op(1,pn)))
    fi
  else eval(f)
  fi
end:

`pd/TD` := proc (f,x,p)
  option remember; 
  if type(x, 'name') then `pd/TD/1`(f,x,p)
  else `pd/TD/1`(evaln(TD(f,`count/r`(x))), `count/f`(x), p)
  fi
end:

`pd/TD/1` := proc (f,x,p)
  option remember; 
  local qs;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
997
  qs := [op(`vars/1`(f))];
Hynek Baran's avatar
Hynek Baran committed
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
  TD(pd(f,p),x) + convert(map(
    proc(q,f,x,p)
      local qp; qp := pd(TD(q,x),p); 
      if qp <> 0 then pd(f,q)*qp else NULL fi end,
    qs, f, x, p), `+`)
end:

#
#   D e r i v a t i v e s   o f   f u n c t i o n s 
#

# For every concrete nary function  a  there should exist a function
# `der/f` of arguments  der,f1,...,fn,x  which returns the derivative
# der of a(f1,...,fn) with respect to x. The actual parameter der will
# eventually be one of pd, TD.

`der/` := proc (der,f,x) der(f,x) end:

`der/Diff` := proc (der,f,x) der(f,x) end:

`der/*` := proc (der,f,x) local i,s,di;
  s := 0;
  for i to nops(f) do di := der(op(i,f),x);
    if di <> 0 then s := s + subsop (i=di, f) fi
  od; s
end:

`der/&*` := proc ()
  local d,f,x,i,s,di;
  d := args[1];
  f := &*(args[2..nargs-1]); x := args[nargs];
  s := 0;
  for i to nops(f) do di := d(op(i,f),x);
    if di <> 0 then s := s + subsop (i=di, f) fi
  od; s
end:

`der/^` := proc (der,f,g,x)
  if type (g,'integer') then g*f^(g-1)*der(f,x)
  else g*f^(g-1)*der(f,x) + ln(f)*f^g*der(g,x)
  fi
end:

`der//unknown` := proc (der,f)  # M.M. 26.11.2004
  local as,x,i,ind;
  global DER;
  if type(f,'indexed') and op(0,f) = DER then 
    ind := op(f);
    if nargs = 4 then
      as := args[3]
    else ERROR(`wrong number of arguments`)
    fi;
    x := args[nargs];
    add(DER[op(sort([i,ind]))](as)*der(op(i,as), x), i = 1..(nargs - 3))
  else
    as := [args[3..(nargs - 1)]];
    x := args[nargs];
    add(DER[i](f(op(as)))*der(as[i], x), i = 1..(nargs - 3))
    fi
end:
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1058
1059
1060
1061
1062
1063
1064
# HB:
# there is an example how to specify derivatives of general function:
`der/MyFun`    := proc (der,f,x) 'D_MyFun'(f)*der(f,x) end:
`der/D_MyFun`  := proc (der,f,x) 'D2_MyFun'(f)*der(f,x) end:
`der/D2_MyFun` := proc (der,f,x) 'D3_MyFun'(f)*der(f,x) end:
`der/D3_MyFun` := proc (der,f,x) 'D4_MyFun'(f)*der(f,x) end:
# :HB
Hynek Baran's avatar
Hynek Baran committed
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092

`der/exp` := proc (der,f,x) exp(f)*der(f,x) end:
`der/ln`  := proc (der,f,x) der(f,x)/f end:

`der/sin` := proc (der,f,x) cos(f)*der(f,x) end:
`der/cos` := proc (der,f,x) -sin(f)*der(f,x) end:
`der/tan` := proc (der,f,x) der(f,x)/cos(f)^2 end:
`der/cot` := proc (der,f,x) -der(f,x)/sin(f)^2 end:
`der/csc` := proc (der,f,x) -cos(f)*der(f,x)/sin(f)^2 end:
`der/sec` := proc (der,f,x) sin(f)*der(f,x)/cos(f)^2 end:

`der/arcsin` := proc (der,f,x) der(f,x)/sqrt(1-f^2) end:
`der/arccos` := proc (der,f,x) -der(f,x)/sqrt(1-f^2) end:
`der/arccot` := proc (der,f,x) -der(f,x)/(1+f^2) end:
`der/arctan` := proc (der,f,x) der(f,x)/(1+f^2) end:

`der/sinh` := proc (der,f,x) cosh(f)*der(f,x) end:
`der/cosh` := proc (der,f,x) sinh(f)*der(f,x) end:
`der/tanh` := proc (der,f,x) der(f,x)/cosh(f)^2 end:
`der/coth` := proc (der,f,x) -der(f,x)/sinh(f)^2 end:
`der/csch` := proc (der,f,x) -cosh(f)*der(f,x)/sinh(f)^2 end:
`der/sech` := proc (der,f,x) -sinh(f)*der(f,x)/cosh(f)^2 end:

`der/arcsinh` := proc (der,f,x) der(f,x)/sqrt(1+f^2) end:
`der/arccosh` := proc (der,f,x) der(f,x)/sqrt(f^2-1) end:
`der/arccoth` := proc (der,f,x) der(f,x)/(f^2-1) end:
`der/arctanh` := proc (der,f,x) der(f,x)/(1-f^2) end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1093
1094
1095
1096
`der/dilog` := proc (der,f,x) # AS
  (ln(f)/(1-f))*der(f,x)
end:

Hynek Baran's avatar
Hynek Baran committed
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
`der/LambertW` := proc (der,f,x)
  LambertW(f)/(1 + LambertW(f))*der(f,x)/f
end:

`der/AiryAi` := proc(der,n,f,x) 
  if nargs = 3 then AiryAi(1,n)*der(n,f) 
  elif n = 1 then AiryAi(f)*f*der(f,x)
  elif pd(n,x) = 0 then AiryAi(n + 1, f)*der(f,x)
  else 'der'(AiryAi(n,f),x)
  fi 
end:

`der/AiryBi` := proc(der,n,f,x) 
  if nargs = 3 then AiryBi(1,n)*der(n,f) 
  elif n = 1 then AiryBi(f)*f*der(f,x)
  elif pd(n,x) = 0 then AiryBi(n + 1, f)*der(f,x)
  else 'der'(AiryBi(n,f),x)
  fi 
end:

`der/RootOf` := proc(der,f,x) 
  if has(f,x) then ERROR(`Sorry, this version of Jets cannot handle`) 
  else 0
  fi 
end:

# Printing partial derivatives

`print/pd` := proc (q,x) 'Diff'(q,`count//seq`(x)) end: 

# Diff is mapped to pd

unprotect(Diff):

Diff := proc(f) 
  if nargs = 1 then
    'Diff(f)' # Maple procedures may use Diff for other purposes
  else pd(f,convert([op(2..nargs, [args])], `*`))
  fi 
end:

#
#   A s s i g n i n g   p a r t i a l   d e r i v a t i v e s 
#

# Values assigned to partial derivatives are stored in a table.

`pd/tab` := table([]):

# Extracting values from the pd table:
# `pd//tab`(f,p) returns the value of pd(f,p).
# If ambiguous, chooses that of minimal size

`pd//tab` := proc (f,p)
  option remember; 
  local el,el1,e,ans,m,q,ql,rt,lb; 
  global `pd/tab`,`dep/tab`; 
  if type (f,'dep') and ars(p) minus `dep/tab`[f] <> {} then RETURN(0) fi;
  lb := `PD:`; rt := `report/tab`[pd];
  if assigned(`pd/tab`[f]) then
    if assigned(`pd/tab`[f][p]) then RETURN(`pd/tab`[f][p]) fi;
    el := op(2,op(`pd/tab`)[f]); # list of already assigned derivatives
    el := map(proc(e,p) p/op(1,e) = op(2,e) end, el, p); # reciprocals
    el := select(proc(e) type (op(1,e),'count') end, el); # select counts 
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1161
    if el = [] then RETURN('pd'(f,p)) fi; # if no counts at all
Hynek Baran's avatar
Hynek Baran committed
1162
1163
    el1 := select(proc(e) not type (op(1,e),{`*`,`^`}) end, el); # 1st order 
    if el1 <> [] then
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1164
      if select(type, el1, anything='numeric') <> [] then ans := 0  
Hynek Baran's avatar
Hynek Baran committed
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
      else
        if rt > 3 then report(lb,[`extending table 1 step: `, nops(el1)]) fi;
        el1 := map(proc(e) Simpl(eval(pd(op(2,e),op(1,e)))) end, el1); 
        if nops(el1) > 1 then
          el1 := sizesort(el1,size);
          ans := op(1,el1); # the smallest size derivative
        else
          ans := el1[1]
        fi
      fi;
      if rt > 3 then report(lb,[`storing a value for`, 'pd'(f,p)]) fi;
      `put/1`('pd'(f,p), ans);
      RETURN (ans)
    else # no 1st order counts
      if select(type, el, anything='numeric') <> [] then
        ans := 0; `put/1`('pd'(f,p), ans); RETURN (ans)
      else
        el := map(proc(e,p) [op(e),`count/length`(op(1,e))] end, el);
        el := sort(el, proc(e1,e2) evalb(op(3,e1) < op(3,e2)) end);
        m := op(3,el[1]); # length of the shortest count
        el := select(proc(e,m) evalb(op(3,e) = m) end, el, m); 
        for e in el do # running through lowest (= m-th) order derivatives
          ql := `pd//tab/1`(op(1,e)); # list of indeterminates
          for q in ql do ans := pd(op(2,e),q);
            if rt > 3 then report(lb,[`storing `, 'pd'(f,p/op(1,e)*q)]) fi;
            `put/1`('pd'(f,p/op(1,e)*q), ans)
          od
        od;
        RETURN(`pd//tab`(f,p))
      fi
    fi 
  fi;
  'pd'(f,p)
end:

`pd//tab/1` := proc(q) # count -> list of indeterminates
  if type (q,`*`) then 
    map(proc(z) if type (z,`^`) then op(1,z) else z fi end, [op(q)]);
  elif type (q,`^`) then [op(1,q)]
  else [q]
  fi
end:

Hynek Baran's avatar
v 5.73a    
Hynek Baran committed
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
`pd/tab/print/unk` := proc(u) # print all assigments of unknown u stored in pd/tab
  global `pd/tab`;
  if assigned(`pd/tab`[u]) then
    map(i -> lprint('pd'(u,op(i))=pd(u,op(i))), [indices(`pd/tab`[u])])
  else
    printf("%a not found in pd/tab\n", u);
  fi
end:


Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1218
1219
1220
1221
1222
1223
1224
1225
#`pd/decomp` := proc(e)
#	 # decomposes pd(f,p) to (f,p) and f to (f,1)
#   if type (e,'name') then (e,1)
#   elif type (e,specfunc(anything,'pd')) then	(op(1,e), op(2,e))
#  else error e, "cannot be decomposed, its not symbol or pd" fi;
#end:


Hynek Baran's avatar
Hynek Baran committed
1226

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1227
1228

# put assigns values to pd's via `pd/tab` or `dep/tab`
Hynek Baran's avatar
Hynek Baran committed
1229
1230
put := proc()
  local e;
Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
1231
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
1232
  if type ([args], 'list'(`=`)) then
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1233
1234
    for e in args do 
    	#`put/items/add`(lhs(e));
Hynek Baran's avatar
v 5.87    
Hynek Baran committed
1235
    	`put/1`(op(e), _rest) 
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1236
    od
Hynek Baran's avatar
Hynek Baran committed
1237
1238
1239
1240
1241
1242
  else ERROR (`wrong arguments`)
  fi;
  `put/evaluate`();
  NULL
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1243

Hynek Baran's avatar
v 5.87    
Hynek Baran committed
1244
unassign('`put/limit/size`', '`put/limit/length`'): # use with extreme care!
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1245

Hynek Baran's avatar
v 5.87    
Hynek Baran committed
1246
1247
1248
1249
1250
1251
`put/1` := proc(p,a) # {`put/limit/size`::{numeric,infinity}:=infinity, `put/limit/length`::{numeric,infinity}:=infinity}
  global `report/tab`, `unk/s`,`unk/<</list`,`dep/tab`,`pd/tab`, `put/name/tab`, `put/limit/size`, `put/limit/length`;
  Report(0, [p]);
  Report(1, [p=size(a)]);
  Report(2, [p=[LVar(a),size(a)]]);
  Report(4, [p=a]);
Hynek Baran's avatar
Hynek Baran committed
1252
1253
1254
  if type (p,'name') then
    `unk/s` := map(`put/1/remove`, `unk/s`, p);
    `unk/<</list` := map(`put/1/remove`, `unk/<</list`, p);
Hynek Baran's avatar
Hynek Baran committed
1255
     `put/name/tab`[p] := a;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1256
    assign(p = a);
Hynek Baran's avatar
Hynek Baran committed
1257
1258
1259
  elif type (p,specfunc(anything,'pd')) then
    if a = 0 and type (op(1,p),'dep') and type(op(2,p),'var') then 
      `dep/tab`[op(1,p)] := `dep/tab`[op(1,p)] minus {op(2,p)}
Hynek Baran's avatar
Hynek Baran committed
1260
    else     	  
Hynek Baran's avatar
v 5.87    
Hynek Baran committed
1261
1262
1263
1264
       if type(`put/limit/length`, numeric) and length(a) > `put/limit/length` then
          printf("put: Ignoring too long (size=%a, length=%a>%a) argument %a=%a\n", size(a), length(a), `put/limit/length`, p, a);   
       elif type(`put/limit/size`, numeric) and size(a) > `put/limit/size` then
          printf("put: Ignoring too large (size=%a>%a, length=%a) argument %a=%a\n", size(a),`put/limit/size`, length(a), p, a);             	  
Hynek Baran's avatar
Hynek Baran committed
1265
       else
Hynek Baran's avatar
v 5.87    
Hynek Baran committed
1266
1267
1268
      	 if type(`put/limit/size`, numeric) or type(`put/limit/length`, numeric) then 
      	   printf("put: Putting (size=%a<=%a, length=%a<=%a) %a=%a\n", size(a),`put/limit/size`, length(a), `put/limit/length`, p, a); 
      	 fi;
Hynek Baran's avatar
Hynek Baran committed
1269
1270
1271
1272
1273
         `pd/tab`[op(1,p)][op(2,p)] := a;
       fi;	
    fi;
  else 
    if simpl(p-a) = 0 then 
Hynek Baran's avatar
Hynek Baran committed
1274
      printf ("ignoring redundant put of %q\n", p);
Hynek Baran's avatar
Hynek Baran committed
1275
    else
Hynek Baran's avatar
v 5.87    
Hynek Baran committed
1276
1277
      lprint(p);
      lprint(a);
Hynek Baran's avatar
Hynek Baran committed
1278
      error "Contradiction in put( %1 = %2 )", p, a;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1279
    fi;
Hynek Baran's avatar
Hynek Baran committed
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
  fi
end:

`put/1/remove` := proc(a,b) if a = b then NULL else a fi end:

`put/evaluate` := proc()
  local e,el,f,ft,q,qs,eq,rpt; 
  global `dep/tab`,`pd/tab`,`nonzero/s`;
# `refresh/1`(`pd//tab`);
  refresh();
  el := op(2,op(`pd/tab`));
  rpt := table([]);
  for e in el do
    f := op(1,e); ft := op(2,e);
    qs := map(proc(x) op(1,x) end, op(2,op(ft)));
    rpt[f] := table([]);
    for q in qs do eq := Simpl(eval(ft[q]));
      if type(q,'var') and eq = 0 and type (f,'dep') then
        `dep/tab`[f] := `dep/tab`[f] minus {q}
      else rpt[f][q] := eq
      fi
    od
  od;
  `pd/tab` := op(rpt);
  nonzero();
  map(
    proc(a)
      if Simpl(eval(a)) = 0 then
        ERROR(`declared nonzero became zero`, a)
      fi
    end, `nonzero/s`);
end:

Hynek Baran's avatar
Hynek Baran committed
1313
1314
1315
1316
1317
1318
1319
1320
1321
`puts/assignments` := proc()
  global `put/name/tab`;
  op(op(op(`put/name/tab`)));  # 'unk_i'=value_i, 'unk_j'=value_j, ... (unordered)
end:

`clear/assignments` := proc({output::identical(expr, eq):=expr})
  description "unassign all symbols assigned by put()";
  global `put/name/tab`;
  local as;
Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
1322
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
1323
1324
  print(`Clearing assignments.`);
  as := `puts/assignments`();
Hynek Baran's avatar
v 5.88    
Hynek Baran committed
1325
1326
  map(unassign@op, [indices(`put/name/tab`)]); # unassign assigned symbols
  `put/name/tab` := table(); # clear the table of assigned symbols
Hynek Baran's avatar
Hynek Baran committed
1327
1328
1329
1330
1331
1332
1333
1334
1335
  if output=eq then 
    return as;
  elif output=expr then
    return op(map(lhs-rhs, [as])) ;
  else 
    error "Unknown output type %1", output; 
  fi;
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360

#`put/items/add` := proc(e)
#  global `put/items/table`;
#  local f, p;
#  f, p := `pd/decomp`(e);
#  `put/items/table`[f] := `put/items/table`[f] union {p};
#end:
#
#`index/put/items/table` := proc(Idx::list,Tbl::table,Entry::list)
#    if (nargs = 2) then # get returns {} if unassigned
#        if assigned(Tbl[op(Idx)]) then Tbl[op(Idx)];
#        else {};
#        end if;
#    else # assign
#        Tbl[op(Idx)] := op(Entry);
#    end if;
#end proc:
#
#`clear/put/items` := proc()
#	 global `put/items/table`;
#	 `put/items/table` := table(`put/items/table`):
#end:
#
#`clear/put/items`():

Hynek Baran's avatar
Hynek Baran committed
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
# To convert the pd table into a list:

pds := proc()
  op(map(`pds/1`, op(2,op(`pd/tab`))))
end:

`pds/1` :=  proc(p)
   op(map(
     proc(e,f) 
       'pd'(f,op(1,e)) = op(2,e) 
     end, op(2,op(2,p)), op(1,p))) 
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1374
# To convert the pd table into a list of expressions while clearing all
Hynek Baran's avatar
Hynek Baran committed
1375
1376
# assignments to pd's:

Hynek Baran's avatar
Hynek Baran committed
1377
1378
`clear/pds` := proc({output::identical(expr, eq):=expr})
  description "unassign all pd's and cc's assigned by put()";
Hynek Baran's avatar
Hynek Baran committed
1379
  local t,aux;
Hynek Baran's avatar
Hynek Baran committed
1380
  global `pd/tab`, `clear/pds/ccclearing`;
Hynek Baran's avatar
v. 5.91    
Hynek Baran committed
1381
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
1382
  print(`Clearing pds.`);
Hynek Baran's avatar
Hynek Baran committed
1383
1384
1385
1386
  reduce();
  t := copy(`pd/tab`);
  `pd/tab` := table([]);
  refresh();
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1387
  #`clear/put/items`():
Hynek Baran's avatar
Hynek Baran committed
1388
1389
1390
1391
1392
1393
1394
1395
1396
  aux := map(`pds/1`, op(2,op(t)));
  aux := sizesort(aux, size);
  if output=eq then 
    return op(aux);
  elif output=expr then
    return op(map(lhs-rhs, aux)) ;
  else 
    error "Unknown output type %1", output; 
  fi;
Hynek Baran's avatar
Hynek Baran committed
1397
1398
1399
1400
end:

reduce := proc()
  local ans1,ans2;
1401
  print(`Reducing...`);
Hynek Baran's avatar
Hynek Baran committed
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
  ans1 := `reduce/pd`();
  ans2 := `unks/assignments`();
  refresh();
  ans1,ans2 
end:

`unks/assignments` := proc()
  global `unk/tab`;
  local aux;
  aux := map(op,[entries(`unk/tab`)]);
Hynek Baran's avatar
Hynek Baran committed
1412
  op(map(proc(f) if f <> eval(f) then f = eval(f) fi end, aux)) # 1 = 'unk_1', 2='unk_2', ... (in the order of unknowns())
Hynek Baran's avatar
Hynek Baran committed
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
end:

`reduce/pd` := proc()
  local e,el,f,ft,p,pl,q,qs,qe,eq,rpt; 
  global `dep/tab`,`pd/tab`;
  el := op(2,op(`pd/tab`));
  rpt := table([]);
  for e in el do
    f := op(1,e); ft := op(2,e);
    pl := map(proc(x) op(1,x) end, op(2,op(ft)));
    if f = eval(f) then
      # remove by independence
      if type(f,'dep') then
        pl := select(proc(p,f) global `dep/tab`;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1427
            evalb(`vars/1`(p) minus `dep/tab`[f] = {}) 
Hynek Baran's avatar
Hynek Baran committed
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
          end, pl, f)
      fi;
      qs := {op(pl)}; qe := {};
      # remove by pd consequence
      for p in pl do
        for q in qs do
          if type (q/p,'count') then qe := qe union {q} fi;
        od;
        qs := qs minus qe
      od
    else 
      # remove by assignment
      qs := {}; eq := eval(f)
    fi;
    # evaluate
    rpt[f] := table([]);
    for q in qs do eq := Simpl(eval(ft[q]));
      if type(q,'var') and eq = 0 and type (f,'dep') then
        `dep/tab`[f] := `dep/tab`[f] minus {q}
      else rpt[f][q] := eq
      fi
    od
  od;
  `pd/tab` := op(rpt);
  pds()
end:

# Checking still unresolved compatibility conditions

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1457
1458
1459
1460
`cc/time` := 0.0 :

# The old cc implementation
# This cc implementation is deprecated but still available, 
Hynek Baran's avatar
v 5.73a    
Hynek Baran committed
1461
1462
# if you need to use it globally, set 
# cc := proc() `cc/old`(args, _options) end:
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1463
1464
# or you may call `cc/old`() at any moment without any side-effects

Hynek Baran's avatar
v 5.73a    
Hynek Baran committed
1465
1466
cc := proc() `cc/new`(args, _options) end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1467
1468
1469
1470
`cc/old` := proc()
  local e,el,f,ft,zl,p,pl,q,i,j,p1,q1,z,ans,sl,aux,rt,lb,time0; 
  global `dep/tab`,`pd/tab`,`cc/s`,`cc/tab`,`cc/time`;
  time0 := time();
Hynek Baran's avatar
Hynek Baran committed
1471
1472
1473
  ans := NULL; 
  zl := [];
  lb := `CC:`; rt := `report/tab`[cc];
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1474
  
Hynek Baran's avatar
Hynek Baran committed
1475
  el := op(2,op(`dep/tab`)); # list of stored dependences
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1476
  if rt > 3 then report(lb,["looking for cc - dep/tab"]) fi;
Hynek Baran's avatar
Hynek Baran committed
1477
1478
1479
1480
1481
1482
  for e in el do 
    f := op(1,e); # an unknown
    ft := op(2,e); # its dependence set
    if rt > 3 then report(lb,[`dependence set:`, 'f' = f, 'ft'= ft]) fi;
    aux := NULL;
    if f <> eval(f) then # f assigned
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1483
      for p in `vars/1`(eval(f)) minus ft do
Hynek Baran's avatar
Hynek Baran committed
1484
1485
1486
1487
1488
1489
1490
        # pd(f,p) should be zero
        aux := aux, [f, eval(f), p, 0, 1, 1] 
      od;
      zl := [op(zl), aux];
      if rt > 2 then report(lb,[f, cat(`ass/dep: `,nops(zl),` c.c.`)]) fi
    fi
  od;  
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1491
1492
  
  if rt > 2 then report(lb,["looking for cc - pd/tab"]) fi;
Hynek Baran's avatar
Hynek Baran committed
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
  el := op(2,op(`pd/tab`)); # list of stored partial derivatives
  for e in el do
    f := op(1,e); ft := op(2,e);
    if rt > 2 then report(lb,['f' = f]) fi;
    pl := map(proc(x) op(1,x) end, op(2,op(ft)));
    aux := NULL;
    if f <> eval(f) then # f assigned
      for p in pl do 
        # pd(f,p) should have the expected value
        aux := aux, [f, f, p, ft[p], 1, `count/length`(p)]; 
      od;
      zl := [op(zl), aux];
      if rt > 2 then report(lb,[f, cat(`ass/pd: `,nops([aux]),` c.c.`)]) fi
    else 
      if type(f,'dep') then # dependence
        for p in pl do
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1509
          if `vars/1`(p) minus `dep/tab`[f] <> {} then
Hynek Baran's avatar
Hynek Baran committed
1510
1511
1512
1513
1514
1515
              # pd(f,p) should be zero 
            if ft[p] <> 0 then
              aux := aux, [f, ft[p], 1, 0, 1, `count/length`(p)]
            fi
          fi; 
          if ft[p] <> 0 then
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1516
            for q in `vars/1`(ft[p]) minus `dep/tab`[f] do
Hynek Baran's avatar
Hynek Baran committed
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
              # pd(f,p*q) should vanish for each q from pd(f,p)
              aux := aux, [f, ft[p], q, 0, 1, `count/length`(p*q)]
            od
          fi
        od
      fi;
      zl := [op(zl), aux];
      if rt > 2 then report(lb,[f, cat(`dep/pd: `,nops([aux]),` c.c.`)]) fi;
      # cross differentiation
      aux := NULL;
      for i to nops(pl) do
        for j to i - 1 do p := op(i,pl); q := op(j,pl);
          if ft[p] <> 0 or ft[q] <> 0 then 
            p1 := `transform/count`(p/q); q1 := `transform/count`(q/p);
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1531
1532
            if p*q1 <> q*p1 then ERROR(`this should not happen`) fi; # NejSpolNas
            aux := aux, [f, ft[p], q1, ft[q], p1, `count/length`(p*q1)] # ii kind
Hynek Baran's avatar
Hynek Baran committed
1533
1534
1535
1536
1537
1538
1539
          fi
        od
      od;
      if rt > 2 then report(lb,[f, cat(`cross-derivatives: `,nops([aux]))]) fi;
      zl := [op(zl), aux];
    fi
  od;  
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1540
  
Hynek Baran's avatar
Hynek Baran committed
1541
1542
  if rt > 1 then report(lb,[cat(`total number: `,nops(zl),` c.c.`)]) fi;
  zl := [op({op(zl)} minus `cc/s`)];
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1543
  inc(`cc/count/total`, nops(zl));
Hynek Baran's avatar
Hynek Baran committed
1544
1545
  if rt > 1 then report(lb,[cat(`of them new: `,nops(zl),` c.c.`)]) fi;
  zl := select(proc(z,m) evalb(op(6,z) = m) end,
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1546
    zl, min(op(map(proc(z) op(6,z) end, zl)))); ## !
Hynek Baran's avatar
Hynek Baran committed
1547
1548
1549
1550
1551
1552
1553
#   zl := select(
#     proc(z,zl) local z1;
#       for z1 in zl do
#         if type (op(6,z1)/op(6,z),'count') then RETURN(true) fi
#       od;
#       false
#     end, zl, zl);
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1554
   inc(`cc/count/computed`, nops(zl));
Hynek Baran's avatar
Hynek Baran committed
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
  if rt > 1 then report(lb,[cat(`of them minimal: `,nops(zl),` c.c.`)]) fi;
  if rt > 1 then report(lb,[`to be computed:`, nops(zl),` c.c.`]) fi; 
  for z in zl do 
    aux := `cc/1`(op(1..5,z)); 
    if aux = 0 then `cc/s` := `cc/s` union {z} 
    else `cc/tab`[f,op(z)] := aux;
      ans := ans, aux
    fi;
  od;
  if rt > 0 then report(lb,[`number of c.c.:`, nops({ans})]) fi;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1565
  inc(`cc/time`, time()-time0);
Hynek Baran's avatar
Hynek Baran committed
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
  {ans}
end:

`cc/s` := {}:

`cc/1` := proc(f,g,p,h,q)
  local ans,pt;
  global `cc/tab`;
  if assigned(`cc/tab`[f][g,p,h,q]) then
    Simpl(eval(`cc/tab`[f][g,p,h,q]));
  elif p = 1 then
    if q = 1 then Simpl(eval(g - h))
    else Simpl(eval(g - pd(h,q)))
    fi
  elif q = 1 then Simpl(eval(pd(g,p) - h))
  else Simpl(eval(pd(g,p) - pd(h,q)))
  fi
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1585
1586
# The common initialization

Hynek Baran's avatar
Hynek Baran committed
1587
1588
`clear/cc` := proc()
  global `cc/s`,`cc/tab`;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1589
  local T;
Hynek Baran's avatar
Hynek Baran committed
1590
  print(`Clearing cc's.`);
Hynek Baran's avatar
Hynek Baran committed
1591
  `cc/s` := {};
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1592
  `cc/tab` := table([]);  
Hynek Baran's avatar
Hynek Baran committed
1593
1594
  CC:-init();
  return NULL;
Hynek Baran's avatar
Hynek Baran committed
1595
1596
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
### New cc iplementation

  CCSidePrice1:= proc(U, p, m)::numeric;
   option remember;
   local q, v, s;
   if `JetMonomTools/J`:-divide1(p,m,q) then
     s, v := `CCSidePrice1/1`(U,p); 
     if nops(`vars/1`(J2j(q)) minus v) > 0 then # trik: pd(U,p) does not depend on q=p/m
       return 0
     else
       return degree(q)*nops(v)*s; 
     fi
   else 
     error "Cannot compute price at %1: %2 |/ %3 ", U, p, m; 
   fi; 
  end:
  
  `CCSidePrice1/1` := proc(U,p)
    option remember;
    local P;
    P := pd(U,J2j(p));
    return (size(P)), `vars/1`(P);
  end:

Hynek Baran's avatar
v 5.73a    
Hynek Baran committed
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701

#`cc/new` := proc({leaveTrivial::truefalse:=evalb(resultType=list), 
#            resultType::identical(set,list):=set, 
#            pop::truefalse:=false,
#            totalNumber::symbol:='None', # output parameter: 
#                              # total number of cc remaining (including the portion returned)
#            ##  keywords for selection of a "portion" of cc's
#            # In all parameters below, -1 means no limitation.
#            # At least 1 cc (if exists) is ALWAYS returned regardless of any selecting criteria.
#            # See CC:-cc for details.
#            maxnum::{posint,identical(-1)} := -1, # total maximum number of returned cc's
#            maxnumP::{posint,identical(-1)} := -1, # return the first maxnumP items and add all cc's of equal prices
#            maxprice::{numeric,infinity} := infinity # maximal price of cc's to be returned
#           }, $)
#  local T, ans, time0, cs, as, i; 
#  global rt, lb, `cc/count/total`, `cc/count/computed`, `cc/time`, cctab;
#  lb := `CC:`; rt := `report/tab`[cc];
#  time0 := time();
#  
#  T := j2J(`cc/pd/listall`()); 
#  
#  forget(CCSidePrice1);
#  cs := CC:-cc(T,  sidePriceFunction=CCSidePrice1,
#	             ':-pop'=pop, 
#	             ':-totalNumber'=totalNumber, # TODO(eff): viz nize
#							 ':-maxnum'=maxnum, ':-maxnumP'=maxnumP, ':-maxprice'=maxprice); 
#  as := map(`cc/ass`, J2j(cs));
#  rprintf(1, ["cc found %a are assembled as %a", cs, as]);
#  as := map(Simpl@eval, as) ; # map(simplify, as, size); ### ???
#  rprintf(1, ["cc found are %a", as]);
#  inc(`cc/count/total`, nops(cs));
#  if leaveTrivial then  
#    ans := as;
#    if rt > 0 then report(lb,[`number of c.c.:`, nops(ans)]) fi;
#  else 
#	  # remove 1=1 unless prohibited by leaveTrivial
#    ans := remove(`=`, as, 0); 
#    if nops(ans)=0 and eval(totalNumber) > nops(cs) then 
#			# at least 1 cc must be returned
#			i := 1;
#			while nops(ans)=0 and eval(totalNumber) > nops(cs) do
#        cs := CC:-cc(T,  sidePriceFunction=CCSidePrice1,
#				             ':-pop'=pop, 
#				             ':-totalNumber'=totalNumber, 
#				             ':-maxnum'=i);		
#				  # TODO(eff): Pravdepodobne bude rychlejsi namisto opakovaneho volani CC:-cc(':-pop'=pop)		
#				  # zavolat jen jednou cc(':-pop'=false), 
#				  # vybrat z vysledku co je potreba  a nakonec zavolat CC:-markFF na zvoleny vyber            
#			  as := map(`cc/ass`, J2j(cs));
#		    as := map(Simpl@eval, as) ; # map(simplify, as, size); ### ???
#			  rprintf(1, ["additional cc found %a and assembled as %a", cs, as]);
#        inc(`cc/count/total`, nops(cs));
#				ans := remove(`=`, as, 0);				
#				if not(pop) then i := i+1 fi;				
#			od;
#			if nops(ans)>0 then
#        if rt > 0 then report(lb,[`nontrivial cc (regardless of selecting criteria given) found`]) fi; 
#      else
#	      if rt > 0 then report(lb,[`None nontrivial cc (regardless of selecting criteria given) found`]) fi; 		  
#      fi;				
#	  else
#      if rt > 0 then report(lb,[`numbers of c.c.: nontrivial`, nops(ans), `trivial`, nops(as)-nops(ans)]) fi;   
#		fi;		
#  fi; 
#  forget(CCSidePrice1);
#  
#  inc(`cc/count/computed`, nops(ans));
#  inc(`cc/time`, time()-time0);
#  rprintf(1, ["There are %a cc at the moment. We have had computed %a of totally %a.", 
#              nops(ans),  `cc/count/computed`, `cc/count/total`]);
#  rprintf(3, ["cc are %q", ans]);
#  
#  if rt > 0 then report(lb,[`total c.c. number: `, `cc/count/total`, 
#                            `computed:`, `cc/count/computed`, `time:`, `cc/time`]) 
#  fi;
#  
#  map(eval, convert(ans,resultType));
#end:


`cc/new` := proc({leaveTrivial::truefalse:=evalb(resultType=list), 
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
            resultType::identical(set,list):=set, 
            pop::truefalse:=false,
            totalNumber::symbol:='None', # output parameter: 
                              # total number of cc remaining (including the portion returned)
            ##  keywords for selection of a "portion" of cc's
            # In all parameters below, -1 means no limitation.
            # At least 1 cc (if exists) is ALWAYS returned regardless of any selecting criteria.
            # See CC:-cc for details.
            maxnum::{posint,identical(-1)} := -1, # total maximum number of returned cc's
            maxnumP::{posint,identical(-1)} := -1, # return the first maxnumP items and add all cc's of equal prices
            maxprice::{numeric,infinity} := infinity # maximal price of cc's to be returned
           }, $)
Hynek Baran's avatar
Hynek Baran committed
1714
1715
  local T, ans, time0, cs, as, i, a, aux, TN, trvcnt, trvcnt0; 
  global rt, lb, `cc/count/total`, `cc/count/computed`, `cc/count/computed/trivial`, `cc/time`, cctab;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1716
1717
1718
  lb := `CC:`; rt := `report/tab`[cc];
  time0 := time();
  
Hynek Baran's avatar
Hynek Baran committed
1719
1720
  if pop<>false then error "Unsupported pop option value %1", pop; fi;

Hynek Baran's avatar
v 5.75    
Hynek Baran committed
1721
1722
  if rt > 1 then
    report(lb, "Dependences and `pd/tab` indices:");
Hynek Baran's avatar
Hynek Baran committed
1723
    try
Hynek Baran's avatar
Hynek Baran committed
1724
1725
      aux := `cc/pd/listall`();
      lprint(map(a -> 'a'(op(vars(a)))=aux[a] , [indices(aux, nolist)]));
Hynek Baran's avatar
Hynek Baran committed
1726
    catch:
Hynek Baran's avatar
Hynek Baran committed
1727
1728
1729
      printf("\n%q\n", StringTools:-FormatMessage( lastexception[2..-1] ));
      lprint("indices(pd/tab)"=indices(`pd/tab`, nolist));
      lprint("`cc/pd/listall`()=", aux);
Hynek Baran's avatar
Hynek Baran committed
1730
    end;
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
1731
1732
  fi;
  
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1733
1734
  T := j2J(`cc/pd/listall`()); 
  
Hynek Baran's avatar
Hynek Baran committed
1735
  if rt > 3 then
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
1736
    report(lb, "Looking for cc for :");
Hynek Baran's avatar
Hynek Baran committed
1737
    map(a->lprint('a'(op(vars(a)))=J2j(T[a])), [indices(T, nolist)]);
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
1738
1739
1740
  fi;
 
  
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1741
  forget(CCSidePrice1);
Hynek Baran's avatar
Hynek Baran committed
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
  
  ### take some portions of cc's until something notrivial is found or no more cc's remains
  ans := [];
  TN := infinity; # last known total number of ccs (unknown yet)
  trvcnt := 0; trvcnt0 := 0; # trivial cc counters (total, last step)
  ###
  ### This is not well tested yet!!! 
  ###
  while nops(ans)=0 and (TN-nops(ans)-trvcnt0)>0 do 
  ###
    cs := CC:-cc(T, 
                 #####sidePriceFunction=CCSidePrice1, #### CCSidePrice1 is not a good idea :(
                 ':-totalNumber'='TN', 
                 ':-maxnum'=maxnum, ':-maxnumP'=maxnumP, ':-maxprice'=maxprice
                 );                
    as := map(`cc/ass`, J2j(cs));
    rprintf(2, ["cc found are assembled as %a", as]);
    as := map(Simpl@eval, as) ; # map(simplify, as, size); ### ???
    rprintf(4, ["cc found %a are assembled and simplifyed as %a", cs, as]);  
    `cc/FF`(cs,as);
    inc(`cc/count/total`, nops(cs));
    
    # remove 1=1 unless prohibited by leaveTrivial
    if leaveTrivial then  
      ans := as;
      trvcnt0 := nops(remove(`=`, as, 0)) - nops(ans);
    else 
     ans := remove(`=`, as, 0); 
     trvcnt0 := nops(as) - nops(ans);	
    fi; 
    inc(trvcnt, trvcnt0);
         
    if rt > 1 then report(lb,[`number of nontrivial cc found: `, nops(ans), 
                          `dropped trivial:`, trvcnt0, 
                          `unprocessed remaining cc`, TN-nops(ans)-trvcnt0]) fi;   
  ###  
  od;
  ###
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1780
1781
  forget(CCSidePrice1);
  
Hynek Baran's avatar
Hynek Baran committed
1782
1783
  ans := map(eval, convert(ans,resultType));
      
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1784
  inc(`cc/count/computed`, nops(ans));
Hynek Baran's avatar
Hynek Baran committed
1785
  inc(`cc/count/computed/trivial`, trvcnt);
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1786
1787
1788
1789
1790
  inc(`cc/time`, time()-time0);
  rprintf(1, ["There are %a cc at the moment. We have had computed %a of totally %a.", 
              nops(ans),  `cc/count/computed`, `cc/count/total`]);
  rprintf(3, ["cc are %q", ans]);
  
Hynek Baran's avatar
Hynek Baran committed
1791
1792
1793
1794
1795
  if rt > 0 then report(lb,[`number of returned c.c. (nontrivial): `, nops(ans), 
                          `of totally remaining cc`, TN-trvcnt0]) fi;             
  rprintf(3, ["cc are %q", ans]);
  if rt > 1 then report(lb,[`total computed c.c. number: `, `cc/count/total`, 
                            `computed in this step:`, `cc/count/computed`, `time:`, `cc/time`]) 
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1796
1797
  fi;
  
Hynek Baran's avatar
Hynek Baran committed
1798
1799
  if totalNumber<>'None' then totalNumber := TN - trvcnt0 fi; # output parameter  
  return map(eval, convert(ans,resultType));
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1800
1801
end:

Hynek Baran's avatar
v 5.73a    
Hynek Baran committed
1802
1803

`cc/FF` := proc(cs, as)
Hynek Baran's avatar
Hynek Baran committed
1804
  description "Mark 0 cc's as fullfilled";
Hynek Baran's avatar
v 5.73a    
Hynek Baran committed
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
  zip(`cc/FF/1`, cs, as)
end:

`cc/FF/1` :=proc(c,a)
  if a=0 then 
    CC:-markFF(c[2],lhs(c[3]),rhs(c[3]),CC:-getHS()[c[1]]);
    return NULL;
  else
    return c
  fi;
end:
    

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
`cc/pd/listall` := proc()
  global `dep/tab`, `pd/tab`;
  local T, e, el, f, ft, ds;
  T := table();
  
  el := op(2,op(`dep/tab`)); # list of stored dependences
  for e in el do 
    f := op(1,e); # an unknown
    ft := op(2,e); # its dependence set
    T[f] := `vars/1`(eval(f)) minus ft;
		if f <> eval(f) then T[f] := T[f] union {1} fi; # f assigned
  od;
  #print('T'=eval(T)); 
  
  el := op(2,op(`pd/tab`)); # list of stored partial derivatives
  for e in el do
    f := op(1,e); # an unknown
    ft := op(2,e); # table of its derivatives
    ds := {op(map(lhs, op(2,op(ft))))};
		#print(DS1, f, ft, ds);
		if type(f,'dep') then
		 ds := map(proc(p) p, 
								       op(`vars/1`(p) minus `dep/tab`[f]),
								       op(`vars/1`(ft[p]) minus `dep/tab`[f])
							 end, ds);
		 #print(DS2, f, ft, ds);
		fi;
    if f <> eval(f) then ds := ds union {1} fi; # f assigned
    if assigned(T[f]) then T[f] := T[f] union ds else  T[f] := ds fi;
  od;
   
  return(eval(T));
end:


`cc/ass` := proc(c::uneval)
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
  global rt, lb;
  local u, r, p, q;
  lb := `CC:`; rt := `report/tab`[cc];
  u := c[1];
  r := c[2];
  p := lhs(c[3]);
  q := rhs(c[3]);
  #printf("cc of %a at %a:\n", u, r);
  #printf("%q\n", [''pd1''( 'PdNE'(u,p), r/p)            =  ''pd1''( 'PdNE'(u,q), r/q)]);
  #printf("%q\n", [''pd1''( `pd/noeval`(u,p), r/p)       =  ''pd1''( `pd/noeval`(u,q), r/q)]);
  #printf("%q\n", [pd1( `pd/noeval`(u,p), r/p)           =  pd1( `pd/noeval`(u,q), r/q)]);
  #printf("%q\n", [`pd/noeval`( `pd/noeval`(u,p), r/p)   =  `pd/noeval`( `pd/noeval`(u,q), r/q)]);

  #pd1( pd1(u,p), r/p) -  pd1( pd1(u,q), r/q);
  #lprint (u,r, p,q);
Hynek Baran's avatar
Hynek Baran committed
1869
  if rt > 3 then 
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
1870
1871
1872
1873
1874
1875
1876
    report(lb, "assembling cc of variable %a(%a) at %a in %a, %a:\n
               pd(%a, %a) = pd( pd(%a, %a), %a) === pd( pd(%a, %a), %a) = pd(%a, %a)", 
                u, vars(u), r, p, q, #[indices(`pd/tab`[u], nolist)],
                [vars(`pd/noeval`(u,p)), Vars(`pd/noeval`(u,p))], r/p, 
                u, p, r/p, 
                u, q, r/q,
                [vars(`pd/noeval`(u,q)), Vars(`pd/noeval`(u,q))], r/q)     
Hynek Baran's avatar
Hynek Baran committed
1877
  elif rt > 2 then 
Hynek Baran's avatar
v 5.75    
Hynek Baran committed
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
    report(lb, "assembling cc of variable %a(%a) at %a in %a, %a:\n
               pd( pd(%a, %a), %a) = pd( pd(%a, %a), %a)", 
                u, vars(u), r, p, q,
                u, p, r/p, 
                u, q, r/q);
   lprint("c = ",c);                
  fi;
  #rprintf(3, ["assembling cc pd( pd(%a, %a), %a) = pd( pd(%a, %a), %a)", u,p, r/p, u,q, r/q]);
  pd1( `pd/noeval`(u,p), r/p) -  pd1( `pd/noeval`(u,q), r/q);  			
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1889
    
Hynek Baran's avatar
v 5.73a    
Hynek Baran committed
1890
1891
    

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
`pd/noeval` := proc(F,U)
  global `pd/tab`;
  local res;
  if U=1 then
    return eval(F)
  else
    res := eval(`pd/tab`[evaln(F)][U], 1); 
    if type(res, indexed) then # not found in `pd/tab`
			pd(F,U);
    else
     res
    fi;
  fi
end:

Hynek Baran's avatar
Hynek Baran committed
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
#
#   S i m p l i f y i n g
#

# Mathematically, f = simpl(f);
# it should be possible to test zero by simpl(f) = 0;
# Finally, simpl(f) = 0 iff numer(simpl(f)) = 0 

unprotect(normal):

simpl := proc(f) normal(f) end:

Simpl := proc(f)
  local V;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1921
1922
1923
1924
1925
1926
1927
1928
  global `Simpl/nocollect/s`;
  if  hasfun(f, TD) then 
    WARNING("Simpl has TD!!!"); 
    simpl(f);
  else
    if nargs > 2 then
      error "Too much arguments"
    elif nargs = 2 then
Hynek Baran's avatar
Hynek Baran committed
1929
1930
      V := args[2];
      assert(type(V,sequential));
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1931
1932
1933
1934
1935
    else
      V := Vars(f) ;
      V := select(proc(v,f) type(f,polynom(anything,v)) end, V, f);
      V := convert(V,set) minus `Simpl/nocollect/s`;
    fi;
Hynek Baran's avatar
Hynek Baran committed
1936
    collect(f, V, distributed, simpl)
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1937
  fi;
Hynek Baran's avatar
Hynek Baran committed
1938
1939
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1940
1941
`Simpl/nocollect/s` := {}:

Hynek Baran's avatar
Hynek Baran committed
1942
1943
1944
1945
#
#   F i n d i n g   d e p e n d e n c e   s e t 
#

Hynek Baran's avatar
v 5.88    
Hynek Baran committed
1946
vars := proc({noWarn::truefalse:=false,  forceError::truefalse:=false})
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1947
1948
1949
  `union`(op(map(`vars/1`, [_rest], _options))); # HB multiple arguments
end:

Hynek Baran's avatar
v 5.88    
Hynek Baran committed
1950
`vars/1` := proc(f, {noWarn::truefalse:=false, forceError::truefalse:=false})
Hynek Baran's avatar
Hynek Baran committed
1951
1952
1953
1954
1955
1956
# option remember;
  if type (f,'constant') then {}
  elif type (f,'name') then 
    if type (f,{`b/var`,`f/var`}) then {f}
    elif type (f,'parameter') then {}
    elif type (f,'dep') then select (type,`dep/tab`[f],'var')   
Hynek Baran's avatar
v 5.88    
Hynek Baran committed
1957
1958
    else if forceError then ERROR ("unknown dependence %1 in vars.", f) fi;
         if not noWarn then WARNING ("unknown dependence %1 in vars.", f) fi; 
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1959
         {f}
Hynek Baran's avatar
Hynek Baran committed
1960
    fi
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1961
  elif type (f,{`+`,`*`,`^`,sequential}) then `union`(op(map(procname,[op(f)], _options))) # HB
Hynek Baran's avatar
Hynek Baran committed
1962
1963
  elif type (f,'function') then 
    if type (f, specfunc(anything,jet)) then {f}
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1964
    elif type (f, specfunc(anything,pd)) then procname(op(1,f), _options)
Hynek Baran's avatar
Hynek Baran committed
1965
    elif type (f, specfunc(anything,TD)) then `vars/TD` (op(f))
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1966
    else `union`(op(map(procname,[op(f)], _options)))
Hynek Baran's avatar
Hynek Baran committed
1967
    fi
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1968
  else error("unexpected object %1 in vars.", f);
Hynek Baran's avatar
Hynek Baran committed
1969
1970
1971
  fi
end:

Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1972

Hynek Baran's avatar
Hynek Baran committed
1973
1974
1975
`vars/TD` := proc (f,x) 
  option remember;
  if type (x,'name') then
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1976
    `vars/1`(f) union `union`(op(map(`vars/TD/1`,`vars/1`(f),x)))
Hynek Baran's avatar
Hynek Baran committed
1977
1978
1979
1980
1981
1982
1983
  else `vars/TD`(f,`count/r`(x)) union
    `union`(op(map(`vars/TD/1`,`vars/TD`(f,`count/r`(x)), `count/f`(x))))
  fi
end:

`vars/TD/1` := proc (q,x) 
  option remember;
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
1984
  `vars/1`(evalTD(TD(q,x))) 
Hynek Baran's avatar
Hynek Baran committed
1985
1986
1987
1988
1989
1990
1991
end:

#
#   J e t s 
#

jet := proc (u::`f/var`,x::`b/count`)
Hynek Baran's avatar
Hynek Baran committed
1992
1993
  option remember, system, threadsafe;
  global `eqn/table`;
Hynek Baran's avatar
Hynek Baran committed
1994
  local y, ans;
Hynek Baran's avatar
v 5.90    
Hynek Baran committed
1995
  if assigned(`eqn/table`[u,x]) then # HB cannot use remembertable to store permanent info
Hynek Baran's avatar
Hynek Baran committed
1996
    RETURN (`eqn/table`[u,x])
Hynek Baran's avatar
Hynek Baran committed
1997
1998
1999
2000
2001
2002
2003
2004
  else 
    for y in `b/var/s` do 
      if type (x/y, 'count') then
        if jet(u,x/y) <> evaln(jet(u,x/y)) then
          RETURN (`simpl/jet`(evalTD(TD(jet(u,x/y),y))))
        fi
      fi 
    od;
Hynek Baran's avatar
v 5.90    
Hynek Baran committed
2005
    RETURN ('jet'(u,sort(x,`b/var/list`)));
Hynek Baran's avatar
Hynek Baran committed
2006
2007
2008
2009
2010
  fi
end:

`simpl/jet` := op(simpl):

Hynek Baran's avatar
v 5.72    
Hynek Baran committed
2011
2012

`jet/gen/all` := proc() # HB
Hynek Baran's avatar
v 5.78    
Hynek Baran committed
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
 description "Generates all internal jets of f up to n-th order";
 global `f/var/list`, `b/var/list`, `eqn/list`;   
 local f, n, cs, es;
 if nargs=1 and type(args[1], {posint,0}) then
   op(map(procname, `f/var/list`, args[1]));
 elif nargs=2 and type(args[1],`f/var`) and type(args[2], {posint,0}) then
   f := args[1];
   n := args[2];
   es := map(e -> [lhs(e)],`eqn/list`);
   es := select(e -> e[1]=f, es);
   es := map(e -> e[2], es); # counts of all assigned jets of f
   cs := [`count/gen/all`(`b/var/list`, n)]; # all counts up ti n-th order
   cs := remove(c-> ormap(e -> divide(c,e), es, c) , cs); # remove assigned jets
   f, op(sort(map2(jet,f,cs), `vars/<<`));
 else
   error "Wrong argument(s) %1", args;
 fi;
Hynek Baran's avatar
v 5.72    
Hynek Baran committed
2030
2031
2032
end:


Hynek Baran's avatar
Hynek Baran committed
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
#
#   E q u a t i o n s
#

equations := proc ()
  global `eqn/list`;
  if not type ([args], 'list'(`=`)) then
    ERROR (`arguments not of type '='`)
  fi;
  `eqn/list` := map (
    proc (e)
      if not type (op(1,e), `j/var`) then 
         ERROR (`not a jet variable on lhs`, op(1,e))
      else (op(op(1,e))) = op(2,e)
      fi
    end, [args]);
  `eqn/list` := 
  map(
    proc(e)
      (op(1,e)) = simpl(eval(op(2,e)))
    end, `eqn/list`);
  refresh (); 
  op(map(proc(e) 'jet'(op(1,e)) = op(2,e) end, `eqn/list`));
end:

`eqn/list` := []:
Hynek Baran's avatar
v 5.90    
Hynek Baran committed
2059
`eqn/table` := table(): # HB
Hynek Baran's avatar
Hynek Baran committed
2060
2061
2062
2063
2064
2065
2066
2067
2068

unprotect(equation):
equation := op(equations):

#
#   R e f r e s h i n g
# 

refresh := proc()
Hynek Baran's avatar
v 5.90    
Hynek Baran committed
2069
  global `eqn/table`, `eqn/list`;
Hynek Baran's avatar
Hynek Baran committed
2070
  `eqn/table` := table():
Hynek Baran's avatar
Hynek Baran committed
2071
  map (`refresh/1`, `refresh/list`);
Hynek Baran's avatar
v 5.90    
Hynek Baran committed
2072
2073
  op(map(proc(e)  assign(`eqn/table`[op(1,e)] = simpl(eval(op(2,e)))); forget(jet); end, # HB prevent remembering of future defined jets
         `eqn/list`));
Hynek Baran's avatar
Hynek Baran committed
2074
2075
end:

Hynek Baran's avatar
v 5.90    
Hynek Baran committed
2076
`refresh/1` := proc(f) forget(f, subfunctions=true, forgetpermanent=true) end:
Hynek Baran's avatar
Hynek Baran committed
2077

Hynek Baran's avatar
v 5.90    
Hynek Baran committed
2078
`refresh/list` := [jet, Simpl, `pd//tab`, Vars, vars, pars, unks, pd]:
Hynek Baran's avatar
Hynek Baran committed
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088

# Symmetries

symmetries := proc ()
  local Q; global `eqn/list`;
  Q := table([args]);
  op (map (proc (e,Q) local lhs,rhs;
    lhs := `dummy/j`(op(1,e)); rhs := op(2,e);
    `symm/1`(lhs,Q) - convert (map (proc (q,rhs,Q) 
      `symm/1`(q,Q) * pd(rhs,q)
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
2089
    end, [op(`vars/1` (rhs) minus `b/var/s`)], rhs, Q), `+`)
Hynek Baran's avatar
Hynek Baran committed
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
  end, `eqn/list`, Q))
end:

`symm/1` := proc (q,Q) 
  if type (q,`f/var`) then Q[q]
  else eval(subsop(0 = TD, 1 = Q[op(1,q)], q))
  fi 
end:

# Universal linearization

lin := proc (f)
  local Q; 
  Q := table([args[2..nargs]]);
  convert (map (proc (q,f,Q) `symm/1`(q,Q)*pd(f,q)
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
2105
    end, [op(`vars/1` (f) minus `b/var/s`)], f, Q), `+`)
Hynek Baran's avatar
Hynek Baran committed
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
end:

# Jacobi bracket

Jacobi := proc (fl,gl)
  local i,n;
  if not type (fl,'list') then ERROR(`invalid argument`, fl) fi;
  if not type (gl,'list') then ERROR(`invalid argument`, gl) fi;
  n := nops(fl);
  if n <> nops(gl) then ERROR(n <> nops(gl)) fi;
  [seq(`Jacobi/1`(fl,gl[i]) - `Jacobi/1`(gl,fl[i]), i = 1..n)]
end:

`Jacobi/1` := proc(fl,g)
  convert(map(
    proc(q,fl,g) 
      local k; global `f/var/list`;
      if type (q,`f/var`) then member(q,`f/var/list`,k); fl[k]*pd(g,q)
      else member (op(1,q),`f/var/list`,k); TD(fl[k],op(2,q))*pd(g,q)
      fi
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
2126
    end, [op(`vars/1`(g) minus `b/var/s`)], fl, g), `+`)
Hynek Baran's avatar
Hynek Baran committed
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
end:

# Conservation laws

laws := proc ()
  local k,neqns,e,lhs,rhs,aux,i,q,P; global `eqn/list`;
  P := table([args]); aux := table ([]); neqns := nops(`eqn/list`);
  for k to `f/dim` do aux[op(k,`f/var/s`)] := 0 od;
  for k to neqns do
    e := `eqn/list`[k]; lhs := `dummy/j`(op(1,e)); rhs := op(2,e);
    aux[op(1,lhs)] := aux[op(1,lhs)]
      + `count/sgn`(op(2,lhs)) * TD(P[k],op(2,lhs));
Hynek Baran's avatar
v 5.7    
Hynek Baran committed
2139
    for q in (`vars/1` (rhs) minus `b/var/s`) do